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The problem of a few interacting fermions in quantum physics has sparked intense interest, par- 
ticularly in recent years owing to connections with the behavior of superconductors, fermionic su- 
perfluids, and finite nuclei. This review addresses recent developments in the theoretical description 
of four fermions having finite-range interactions, stressing insights that have emerged from a hyper- 
spherical coordinate perspective. The subject is complicated, so we have included many detailed 
formulas that will hopefully make these methods accessible to others interested in using them. The 
universality regime, where the dominant length scale in the problem is the two-body scattering 
length, is particularly stressed, including its implications for the famous BCS-BEC crossover prob- 
lem. Derivations and relevant formulas are also included for the calculation of challenging few-body 
processes such as recombination. 
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I. INTRODUCTION 



The problem of four interacting particles in nonrelativistic quantum mechanics arises in a number of different 
physical and chemical contexts. P~H1] While tremendous theoretical progress has been achieved in the three-body 
problem, [5 15 particularly in the past two decades, the four-body problem remains still in its infancy by comparison. 
Like the three-body problem, the four-body problem consists of two qualitatively different subcategories, one in which 
some of the particles have Coulombic interactions. [H)H2"3"] and the other subcategory in which all forces between 
particles have a finite range or else have at most a rapidly-decaying multipole interaction at long-range. [TJ [H l2"4T - f2"5] 
The subject of the present review concerns the latter category, which is particularly relevant to modern day studies 
of ultracold quantum gases composed of neutral atoms and/or molecules. The scope of this subject is much broader 
than that of ultracold gases alone, however, as 4-body reactive processes such as AB+CD— s-AC+BD, or — s-A+BCD, 
or — s-A+B+C+D occur in nuclear and high-energy physics as well as in chemical physics. The time-reverse of these 
processes is also important for understanding the loss rate in a degenerate quantum gas, notably the process of 
four-body recombination which had hardly received any attention until very recently. 

While of course many important advances have been achieved in few-body physics without the use of hyperspherical 
coordinates, treatments using these coordinates have real advantages for a number of problems. Early on, for instance, 
Thomas |30j proved an important theorem about the nonzero range of nucleon-nucleon forces, using an analysis in 
which the hyperradial coordinate played a crucial role although he did not refer to it by that name. (See, for instance, 
Eq.lllc of jj.) Further developments in the use of hyperspherical coordinates in collision problems were pioneered 
by Delves [3TJ El] and they played a key role in the derivation of the Efimov effect [33J El] As we will see below, the 
advantages accrue not only in terms of computational efficiency, but also in terms of the insights and quasi-analytical 
formulas that can be deduced for scattering, bound, and resonance properties of the system. For this reason, the 
present review concentrates on the hyperspherical studies of the four-body problem, concentrating on recent progress 
and results that have emerged, and on problems that currently seem ripe for pursuit in the near future. 

In early studies [35, 36 , hyperspherical coordinates were viewed as capable of providing a deeper understanding 
of the nature of exact bound state solutions, for instance for the helium atom [37] ■ And Delves [3T1 132"] used these 
coordinates to discuss rearrangement nuclear collisions from a formal perspective. But a turning point in the utility 
of hyperspherical coordinate methods was introduced by Macek in 1968[3H], in the form of two related tools: the 
adiabatic hyperspherical approximation and the (in principle exact) adiabatic hyperspherical representation. Both 
of these methods single out a single collective coordinate for special treatment, the hyperradius R of the TV-body 
system, which is handled differently from all remaining space and spin coordinates, fi. The hyperradius is a positive 
"overall size coordinate" of the system, whose square is proportional to the total moment of inertia of the system, 
i.e. R 2 = -ryEjmjr?, where rrii is the mass of the i-th particle at a distance rj from the center of mass, and M is any 
characteristic mass which can be chosen with some arbitrariness. |39j . 
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In Macek's adiabatic approximation, the Hamiltonian is diagonalized at fixed values of R, and the resulting energies 
plotted as functions of the hyperradius can be viewed as adiabatic potential curves U^{R) as in the ordinary Born- 
Oppenheimer approximation for diatomic molecules. The first prominent success of the adiabatic approximation was 
the grouping together of He autoionizing levels having similar character into one such potential curve. [38] Subsequent 
studies showed that He and H photoabsorption is dominated by a small subset of such potential curves. [4T)1 - |4"2"] 
suggesting that Macek's adiabatic scheme is much more than just a mathematical technique for solving the Schrdinger 
equation, but that it also provides an insightful physical and intuitive formulation that can be used qualitatively and 
semiquantitatively in the same manner as the Born-Oppenheimer treatment which has been so successful in molecular 
physics. 

At the same time, however, subsequent applications of the strict adiabatic hyperspherical approximation showed 
its limitations. [131 23] Some classes of energy levels or low-energy scattering properties could be described to semi- 
quantitative accuracy, but in other cases it failed to give a reasonable description of the spectrum, sometimes even 
qualitatively. As this has become more and more appreciated, it has become increasingly common to treat few- 
body systems using the adiabatic hyperspherical representation, in principle an exact theory that does not make 
the adiabatic approximation; in this method several adiabatic hyperspherical states are coupled together and their 
nonadiabatic interactions are treated explicitly. Implementation of the adiabatic hyperspherical representation is 
sometimes carried out in exact numerical calculations [8, 45, 46 , but in many cases semiclassical theories such as the 
Landau-Zener-Stueckelberg formulation are sufficiently accurate and useful. |47 j 

In the four-body problem, some initial studies using hyperspherical coordinates were carried out for the description of 
3-electron atoms such as Li, He - , and H -0HIH2] But the method was improved to the point of being a comprehensive 
approach by Refs. [16U22] Despite our focus in the present review article on four interacting particles with short-range 
interactions, we summarize briefly the headway that has previously been achieved for Coulombic systems. For three- 
electron atoms, the topology is of course quite different and more interesting than for two-electron atoms. For 
instance, whereas one observes one or more two-electron hyperspherical potential curves that converge at R — ¥ oo to 
every possible one-electron bound state, the three-electron atom potential curves converge also to unstable resonance 
levels of the residual two-electron ion that have a nonzero autoionizing decay width. There are multiple families of 
potential curves that represent new physical processes such as post-collision interaction in addition to the triply-excited 
states and their decay pathways. Tremendous technical challenges were overcome in an impressive series of articles 
by Lin, Bao, Morishita, and their collaborators, to enable the calculation of accurate hyperspherical potential curves 
for three-electron atoms. [I6-22] For a recent broader review of triply-excited states that also discusses alternative 
approaches beyond the hyperspherical analysis, see [2"3"] . 

Another theoretically challenging type of four-body problem in chemical physics has been the dissociative recombi- 
nation of Hj^ induced by low energy electron collision. Here the 3 bodies are the nuclei (augmented by two "frozen" 
Is electrons that play no dynamical role at low energies), while the fourth body is the incident colliding electron. 
The solution of this problem, including the identification of Jahn- Teller coupling as the controlling mechanism, has 
been greatly aided by the use of hyperspherical internuclear coordinates. They allowed a mapping of the dynamics to 
a single hyperradius, in addition to multichannel Rydberg electron dynamics that could be efficiently handled using 
multichannel quantum defect techniques and a rovibrational frame transformation. [50TI52) 

More relevant to the present review of four-body interactions of short-range character are some long-standing 
problems of reactive processes in nuclear physics and in chemical physics. Fundamental groundwork was laid by 
Kuppermann [551 154] and by Aquilanti and Cavalli[55|. which concentrated on developing coordinate systems and useful 
solutions of the noninteracting problem, which are the hyperspherical harmonics. However, whereas hyperspherical 
harmonics constitute a complete, orthonormal basis set in general, which have numerous useful formal properties, in 
our experience they provide poor convergence when used alone as a basis set to expand a reactive collision wavefunction. 

The tremendous growth of ultracold atomic physics has stimulated much of the current interest in few-body and 
many-body processes that are deeply quantum mechanical in nature. And indeed, some of the progress can be traced to 
the advances that have been made in our understanding of few-body collisions and resonances in the low-temperature 
limit. Some of the most important advances were the development of accurate theoretical models for atom-atom 
collisions at sub-millikelvin temperatures. [55M62] Ab initio theory was not sufficiently advanced to predict the atom- 
atom interaction potentials to sufficient accuracy, so refinements and adjustments of a small number of parameters (the 
singlet and triplet scattering lengths and in some cases the van der Waals coefficient and the total numbers of singlet 
and triplet bound levels) were needed to specify the two-body models. Once the two-body interactions were well 
understood, the next challenge became three-body collisions. In most degenerate quantum gases created during the 
past decade or longer, the lifetimes have been controlled by three-body recombination, i.e. in a single-component BEC, 
this is the process A + A + A — y A%+ A. Advances in understanding and in the ability to carry out nonperturbative 
three-body recombination calculations resulted, by the late 1990s, in some of the first survey studies of the dependence 
of the three-body recombination rate K 3 on the two-body scattering length. Two independent treatments utilizing 
the adiabatic hyperspherical representation[Tni E] led to the prediction that destructive interference minima should 
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exist at positive atom-atom scattering lengths a, with universal scaling behavior connected intimately with the Ehmov 
effect. Such minima have apparently been observed recently in experiments. 63J Ref.^lO additionally predicted that 
three-body shape resonances, also connected intimately with the Ehmov effect, should arise periodically in a, and the 
first such Efimov resonance was observed experimentally in 2006 by the Innsbruck group of Grimm. [64] 

Not long after the dependence of K3 on a had been identified by the aforementioned theoretical treatments in hyper- 
spherical coordinates, alternative treatments provided different ways to understand many of these results. Effective 
field theory [63], functional renormalization|66] . Faddeev treatments in momentum space[67 ( 68 , a transition matrix 
approach based on the three-body Green's function [SS], an d an analytically-solvable model treatment of the Efimov 
problem [70] This large number of independent theoretical formulations, which by and large reproduce and in some 
cases extend the 1999 predictions, is an encouraging confluence that suggests our understanding of the three-body 
problem with short-range forces is nicely on track. 

In contrast, the description of many four-body scattering processes, especially those with a final or initial state 
having four free particles such as the recombination process A + A + A + A — > A3 + A or A2 + A2 or A2 +A + A, is a field 
in its infancy by comparison with the state of the art for the three-body problem. Most previous attention to date has 
concentrated on either four-body bound states such as the alpha particle ground or excited states [7Tj. or else simple 
exchange reactions with two-body entrance and exit channels, such as H + H.%0 — s- H2 + OH [HH]. Some theoretical 
results of this class have been derived in the context of ultracold fermi or bose gases. One of the most important was 
the prediction by Petrov, Salomon, and Shlvapnikov|72j that the rate of inelastic collisions, between weakly bound 
dimers composed of two equal mass but opposite spin fermions, should decay at large fermion-fermion scattering 
lengths as a~ 2 - 55 . Two experiments are consistent with this prediction. [731 HI] This result has been confirmed at 



a qualitative level in a separate hyperspherical coordinate treatment discussed below in Sec. VD but with some 
quantitative, temperature-dependent differences. [75] The real part of the scattering length, associated with purely 
elastic scattering, is also important in the BEC-BCS crossover problem, and its value has been predicted in a number 
of independent studies to equal 0.6a. 

For four identical bosons with large scattering lengths, an insightful theoretical conjecture by Hammer and 
Platter[75] suggested that two four-body bound levels should exist that are attached to every 3-body Efimov state. 
When this problem was tackled using the toolkit of hyperspherical coordinates, the resulting potential curves and 
their bound and quasi-bound levels provided strong numerical evidence in support of this conjecture. [77] Moreover, 
once the major technical challenge of computing the adiabatic hyperspherical potential curves had been overcome, 
through the use of a correlated Gaussian basis set expansion [75], it was possible to calculate four-body recombination 
rates and demonstrate that signatures of four-body physics had in fact already been present and observed in the 2006 
Efimov paper by the experimental Innsbruck group.[64 The four-body resonance features had not been interpreted 
as such in that study, but a subsequent experiment by the same group [79] provided strong confirmation of this point. 
This theoretical development was also aided by a general derivation of the N-body recombination rate [50] in terms of a 
scattering matrix determined within the adiabatic hyperspherical representation. Further extensions have permitted 
an understanding of dimer-dimer collisions involving four bosonic atoms. [81] Another positive advance during the 
last few years has been a treatment of atom-trimer scattering that has determined the lifetime of universal bosonic 
tetramer states [82] and the analysis of the Efimov trimer formation via four-body recombination. [83] 

Our aims in this review are to present some of the technical developments that have recently enabled an extension 
of the adiabatic hyperspherical framework that can handle four or more particles. The most technically-challenging 
aspect of this is the solution of the fixed-hyperradius Schrodinger equation to determine the adiabatic potential 
curves and their couplings that drive inelastic, nonadiabatic processes. Once those couplings and potential curves 
are known, it is comparatively simple and intuitive to understand at a glance the competing reaction pathways that 
can contribute to any given process. In many cases, those pathways are sufficiently small in number, and sufficiently 
localized in the hyperradius, to permit semiclassical WKB and Landau-Zener-Stueckelberg-type theories[I71 154T456] 
to give a semiquantitative description. Such approximate treatments are especially useful for interpreting the results 
of quantitatively accurate coupled channel solutions to the coupled equations. 



II. GENERAL FORM OF THE ADIABATIC HYPERSPHERICAL REPRESENTATION 



One of the greatest advantages of using the hyperspherical adiabatic representation is that it offers a simple, 
yet quantitative, picture of the bound and quasi-bound spectrum as well as scattering processes. It reduces the 
problem to the study of the hyperradial collective motion of the few-body system in terms of effective potentials 
and where inelastic transitions are driven by non-adiabatic couplings. The effective potentials, and the couplings 
between different channels, offer a unified, conceptually clear, picture of all properties of the system. Below, we give 
a general description of the adiabatic hyperspherical representation for a general iV-body problem. Details regarding 
the Coordinate transformations that accomplish the conversion from Cartesian to angular variables (along with R) 
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are given in Appendix [A"| 

A. Channel Functions and Effective Adiabatic Potentials 



In the adiabatic hyperspherical representation, the iV-body Schrodinger equation can be written in terms of the 
rescaled wave function ip = i?( 3 ( JV ~ 1 )~ 1 )/ 2 v[f (hi atomic units), 



1 d 2 
2/idR 2 



ip(R,n) = Etp(R,n), 



(i) 



where \i is the arbitrary, reduced, mass and E is the total energy. It is interesting to notice that the above form of 
the Schrodinger equation is the same irrespective of the system in question, leaving all the details of the interactions 
in the adiabatic Hamiltonian H ad (R, f2), where VI denotes the set of all hyperangles. 

The few-body effective potentials are eigenvalues of the adiabatic Hamiltonian , obtained for fixed values of R, 
i.e., with all radial derivatives omitted from the operator: 



H ad (R,n)^{R;n) = U u (R)<f> u (R; fi). 



(2) 



where & U (R; f2), the eigenstates, are the channel functions, U V (R) the few-body potentials, and the adiabatic Hamil- 
tonian given by 



H ad {R,n) 



A 2 (n) + (3N - 4)(3iV - 6)/4 
2fiR 2 



v(R, n). 



(3) 



The operator A is the squared grand angular momentum defined in Eq. ( B3 1 , and V contains all the interparticle 
interactions. In the above equations, v is a collective index that represents all quantum numbers necessary to label 
each channel. 

From the analysis above, since the channel functions ^ V {R] f2) form a complete set of orthonormal functions at each 
R, they are a natural base to expand the total rescaled wavefunction, 



(4) 



where the expansion coefficient F V (R) is the hyperradial wave function. In this representation, the total wave function 
is, in principle, exact. Upon substituting Eq. Q into the Schrodinger equation ([TJ and projecting out <&„, the 
hyperradial motion is described by a system of coupled ordinary differential equations 



1 d 2 
2(1 dR 2 



U V {R) 



F„{R) 



-E 



2P VV ,(R)—+Q VV ,{R) 
dR 



F V ,(R)=EF V {R\ 



(5) 



where P vv i (R) and Q vv i (R) are the nonadiabatic coupling terms responsible for the inelastic transitions in A^-body 
scattering processes. They are defined as 



and 



/',„,(./.') = ((«!>,.( /?.<2) ^ Sv(/M>! 



(6) 



Qvu>{R) 



dR 2 



(7) 



where the double brackets denote integration over the angular coordinates £1 only. 

Although in the adiabatic hyperspherical representation the major effort is usually in solving the adiabatic equation 
([2]), the hyperradial Schrodinger equation ^ is central to the simplicity of this representation. Since R represents 
the overall size of the system, the hyperradial equation ^ describes the collective radial motion under the influence 
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of the effective potentials W v , denned by 



W V (R) = U V {R) 



1 



Qw{R), 



(8) 



2/t 



while the inelastic transitions are driven by the nonadiabatic couplings P vv i and Q vv ' . Scattering observables, as well 
as bound and quasi-bound spectrum, can then be extracted by solving Eq. ([5]). As it stands, Eq. ^ is exact. In 
practice, of course, the sum over channels must be truncated, and the accuracy of the solutions can be monitored 
with successively larger truncations. Therefore, in the adiabatic hyperspherical representation the usual complexity 
due to the large number of degrees of freedom for few-body systems is conveniently described by a one-dimensional 
radial Shcrodinger equation, reducing the problem to a "standard" multichannel process. 

The hyperspherical adiabatic representation has been shown to offer a simple and unifying picture for describing 
few-body ultracold collisions in the regime where the short-range two-body interactions are strongly modified due 
to a presence of a Fano-Feshbach resonance [57]. In this regime, the long-range properties of the few-body effective 
potentials W v become very important and other analytical few-body collision properties can be derived. For instance, 
the asymptotic behavior of the few-body effective potentials W v determine the generalized Wigner threshold laws 
for few-body collisions [88] , i.e., the energy dependence of the ultracold collisions rates in the near-threshold limit. 
Moreover, when the two-body interactions are resonant, few-body effective potentials are modified accordingly to 
universal physics [BUI HO], as we will show in the following sections. From this analysis, a simple picture describing 
both elastic and inelastic transitions emerges. We also discuss the validity of our results in the context of numerical 
calculations carried out through the solution of Eqs.pl and (p>]) for a model two-body interaction. 



Here we derive a formula for the generalized cross-section describing the scattering of N-particlcs. Our formulation 
is based on the solutions of the hyperadial equation ^ but is sufficiently general to describe any scattering process 
with particles of either permutation symmetry. The only information required in our derivation is that at large R, 
the solutions to the angular portion of the Schrodinger equation yield the fragmentation channels of the TV-body 
system, i.e., the same asyptotic form of the adiabatic effective potentials p]), and the quantum numbers labeling those 
solutions index the 5-matrix. 

This derivation begins by considering the scattering by a purely hyperradial finite-range potential in <i-dimensions, 
then the resulting cross section is generalized to the case of anisotropic finite range potentials in (i-dimensions "by 
inspection " , which we interpret in the adiabatic hyperspherical picture. For clarity, we adopt a notation that resembles 
the usual derivation in three dimensions. 

In (i-dimensions, the wavefunction at large R behaves as: 



B. Generalized Cross Sections 



* 7 -> e 



ik-R 



+ f(k,k') 



R(d-l)/2 



(9) 



Equivalently, an expansion in hyperspherical harmonics is written in terms of unknown coefficients Ax^: 




(10) 



Here, Y\u are hyperspherical harmonics and jf (nf) are hyperspherical Bessel (Neumann) functions [91] . 




(11) 



where a = d/2 — 1. We will make use of the asymptotic expansion, 




(12) 
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and the plane wave expansion in d-dimensions |91j : 



ik-R 



( d 2)!!— — ^ Z A jt(^)^(fc)n,( J R). (13) 



r(d/2) A ; , 



Identifying the incoming wave parts of 1 $> I and ty 11 yields the coefficients A^^: 

A Xv = e iS ^d-2)\\^^i x Y^(k). (14) 

Inserting the coefficients A\ v back into the expression for ty 11 gives the expression for the scattering amplitude: 

d-i 

f(k 1 k')=f 2 ^) 2 Y,n,(k)Y x ,(k')(e^-l). (15) 
The immediate generalization of this elastic scattering amplitude to an anisotropic short-range potential is of course: 

d-l 

2tt n 



/(fe,fe') = ( — ) n;(fe)n v (fc')(^,A>' - 5 AV <W). (16) 

A/iA' fi' 



Upon integrating \f(k, k')\ 2 over all final hyperangles k, and averaging over all initial hyperangles k as would be 
appropriate to a gas phase experiment, we obtain the average integrated elastic scattering cross section by a short- 
range potential: 



7 v ' A/xA'/x' 



where f2(d) = 2ir d / 2 /T(d/2) is the total solid angle in <i-dimensions $T\. This last expression is immediately interpreted 
as the average generalized cross section resulting from a scattering event that takes an initial channel into a final 
channel, i = X'fjf — > X/i = f. Since this S-matrix is manifestly unitary in this representation, it immediately applies 
to inelastic collisions as well, including iV-body recombination, in the form: 

'■"Hir^.-*/.! 5 . <i8) 

It is worth noting that this expression needs to be simply summed up for all initial and final channels contributing 
to a given process of interest, including degeneracies. For instance, we note that in the case of a purely hyperradial 

potential, each A has Mid. A) = - - ^ , \ — i^-h — 1 degenerate values of u. 

P v ' ' r(A + l)r(d- 1) B P 

In this form, we can readily interpret the generalized cross section derived above in terms of the unitary S-matrix 
computed by solving the exact coupled-channels reformulation of the few-body problem in the adiabatic hyperspherical 
representation 38 . In principle this can describe collisions of an arbitrary number of particles. Identical particle 
symmetry is handled by summing over all indistinguishable amplitudes before taking the square, averaging over solid 
angle, then integrating over distinguishable final states to obtain the total cross section: 

indist f dk r dk « .-,, |2 , r dist 

a = jN p J W) 1 Pfi ' )l P ■ 

Here N p is the number of terms in the permutation symmetry projection operator (e.g. for N identical particles, 
N P = M.) 

The cross section for total angular momentum J and parity II includes an explicit 2 J + 1 factor. Hence, the cross 
section from incoming channel i to the final state /, properly normalized for identical particle symmetry, is given in 
terms of general S'-matrix elements as [80] 



idist ( ril 



27T\ d - 1 1 



(■n = N p[j -) pj + DSfr-Sfi'. (19) 
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For the process of N-body recombination in an ultracold trapped gas that is not quantum degenerate, the experi- 
mental quantity of interest is the recombination event rate constant Kff which determines the rate at which atoms 
are ejected from the trapping potential: 

N=2 v ; 

where n is the number density. The above relation assumes that the energy released in the recombination process is 
sufficient to eject all collision partners from the trap. The event rate constant [recombination prob abil ity per second 



for each distinguishable TV-group within a (unit volume)^ ^] is the generalized cross section Eq. (191 multiplied by 



a factor of the iV-body hyperradial "velocity" (including factors of h to explicitly show the units of Kn): 

k£ =J2—*7l dis V u )- (2i) 

Here, the sum is over all initial and final channels t hat contribute to atom-loss. 

The relevant S-matrix element appearing in Eq. ( 19 ) from an adiabatic hyperspherical viewpoint is Sjj where i 



and / are the initial and final channels (i.e. solutions to Eq. ([2| in the limit R — > oo). In the ultracold limit, the 
energy dependence of the recombination process is controlled by the long-range potential Eq. Q in the entrance 
channel i — > A m j ra , where A m j„ is the lowest hyperangular momentum quantum number allowed by the permutation 
symmetry of the iV-particle system. For any combination of bosons and distinguishible particles, A m ;„ = 0, while for 
fermions, the permutation antisymmetry adds nodes to to the hyperangular wavefunction leading to \ m in > 0. 
As a concrete example, consider the recombination formula for the four-fermion process F+F' +F+F' — » FF' +FF' . 



In applying the permutation symmetry operator, it is convenient to employ the H-type coordinates given in Eq. ( C3 1 



Expressing the hyperangular momentum operator in these coordinates, it is possible to show (see Section III B ) that 



A = (I i + I2 + 2ni ) + I3 + 2n2, where l\, I2 and I3 are the angular momentum quantum numbers associated with the 



three Jacobi vectors in Eq. ( C3 ) and ni, n 2 are both non- negative integers. Antisymmetry under exchange of identical 
particles in these coordinates implies that l\ and I2 must be odd. The lowest allowed values are then l± = I2 = 1, 
I3 = rii = 112 = such that X m in = 2. 

The preceding arguments enable us to calculate the generalized Wigner threshold law for strictly four-body re- 
combination processes where the four particles undergo an inelastic transition at R ~ a; any nonadiabatic couplings 
Eqs. ^ and ^ at R ^> \a\ can be viewed as three-body processes with the fourth particle acting as a spectator. 
The asymptotic (R 3> |o|) form of the effective potential can be written in terms of an effective angular momentum 
quantum number l e 

W ^ W ^ with l e = (2X mm + d - 3)/2. (22) 

It was shown in 80J following the treatement of Berry [52] , the WKB tunneling integral gives the threshold behavior 
of the 5-matrix element Sfi cx e^~ 27 ^ with 

/•(3iV-5+2A m j„)/2fc 

7 = Im / dR y/2fi N (E-W'{R)) (23) 

1/4 



and where W'(R) = W{R) + 2i i n r? IS t ne effective potential with the Langer correction [93]. The lower limit of 
the integral R* coincides with the maximum of the nonadiabatic coupling strength Pj i /\Ui(R) — Uf(R)\ defined in 
Eq. ([6]). For recombination into weakly bound dimers or trimers (of size |o|), R* ~ \a\ so that in the threshold limit 

[SO] 

1 - l^l 2 oc e-^ = (k\a\) 2X — +3N - 5 . (24) 
Unitarity of the S'-matrix implies that 1 — \Su\ 2 = J2f^i \Sfi\ 2 i which is related to the total inelastic cross section 



through Eq. (19). If inelastic transitions are dominated by recombination then the scaling law for the recombination 



event rate constant is: 

^cxfc 2A """|a| 2A — +3JV - 5 (25) 
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a-dependence of K± 
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+ B 
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BBB 
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constant 
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F + 


F' + F' -)• 


FF 1 + 


FF' 


2 


E 2 


H 11 


F - 


f F - 


fi + y^ 


FFX 


+ Y 


1 


E 





TABLE I: The generalized Wigner threshold laws are given for a limited set of four-body recombination processes. Here, B 
denotes a boson, F and F' are fermions in different "spin" states, and X and Y are distinguishible atoms. Note that since in 
general the scattering lengths for the F — X , F — Y and X — Y interactions are different, the scaling with respect to a is not 
given for this 3-component case. 

We stress that the above expression gives the overall scaling of the event rate, and in cases where the coupling to 
lower channels occurs in the region r -C R* <C \a\, one must include the additional WKB phase leading to a modified 
scaling with respect to a. The k dependence arises through the outer turning point limit in the WKB tunneling 
phase integral. This occurs, for example, in the case of four- identical bosons treated in [80]. For the four-fermion 
problem, the effective angular momentum quantum numbers for the universal potentials in the region r$ -C R -C \a\ 
are calculated in Ref. |M|- Table | gives the value of \ m in along with the overall recombination rate scaling with \a\ 
for a few select cases. 



III. VARIATIONAL BASIS METHODS FOR THE FOUR-FERMION PROBLEM 

Solution of the four-body hyperangular equation, Eq. poses significant challenges, since the difficulty grows 
exponentially with the number of particles. For four particles with zero total angular momentum, Eq. ^ consist of 
a 5-dimensional partial differential equation. Some state-of-the-art methods for three-particle systems often employ 
B-splines or finite elements. In fact, if 40 — 100 B-splines were used in each dimension to solve Eq. ^ (a common 
number in three-body calculations [THl EH there would be 10 8 — 10 9 basis functions resulting in 10 11 — 10 13 

non-zero matrix elements in a sparse matrix. The computational power required for such a calculation is currently 
beyond reach. Therefore, in order to proceed numerically, a different strategy must be developed. In this review, we 
describe two of our current numerical techniques. 

The method of this section for a two-component system of four fermions uses a non-orthogonal variational basis 
set consisting of some basis functions that accurately describe the system at very large hyperradii, R S> |a|, and 
other functions that describe the system at very small hyperradii, R <C \a\. If both possible limiting behaviors are 
accurately described within the basis, then a linear combination of these two behaviors might be expected to describe 
the intermediate behavior of the system. |42l !Ki 

As with the correlated Gaussian method of Section |IV A| the use of different Jacobi coordinates plays a central role 
in the variational basis method. Depending on the symmetries, interactions, and fragmentation channels inherent in 
the problem, different coordinates may significantly affect the ease with which the problem can be described. For 
example, in the four fermion problem, the fermionic symmetry of the system can be used to significantly reduce the 
size of the basis set needed to describe the possible scattering processes. Describing this symmetry in a poorly-chosen 
coordinate system can create considerable difficulty. The two main types of Jacobi coordinate systems are called H- 
type and K-type, shown schematically in Fig. [I] We discuss some of the relevant properties of the different coordinate 
systems here. Appendix [C] gives a detailed account of the Jacobi coordinate systems used in this review and of the 
transformations between them. 

j k j k 




i I i I 

H-type K-type 



FIG. 1: The two possible configurations Jacobi coordinates in the four-body problem are shown schematically. 
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H-type Jacobi coordinates are constructed by considering the separation vector for a pair of two-body subsystems, 
and the separation vector between the centers of mass of those two subsystems. Physically, H-type coordinates are 
useful for describing correlations between two particles, for example a two-body bound state or a symmetry between 
two particles, or two separate two-body correlations. K-type Jacobi coordinates are constructed in an iterative way 
by first constructing a three-body coordinate set as in Eq. (C2|, and then taking the separation vector between the 
fourth particle and the center of mass of the three particle sub-system. When two particles coalesce (e.g. when r,; = rj 
in Eq. CI), the H-type coordinate system reduces to a three-body system with two of the four particles acting like a 
single particle with the combined mass of its constituents. Locating these "coalescence points" on the surface of the 
hypersphere is crucial for an accurate description of the interactions between particles, and this coordinate reduction 
will prove useful for the construction of a variational basis set. 

Examination of Fig. [T] shows that K-type Jacobi coordinate systems are useful for describing correlations between 
three particles within the four-particle system. In the four-fermion system, there are no weakly-bound trimer states, 
whereby K-type Jacobi coordinates will not be used here, but the methods described in this section can be readily 
generalized to include such states. Unless explicitly stated, all Jacobi coordinates from here on will be of the H-type. 

The task of parametrizing the 3 Jacobi vectors in hyperspherical coordinates remains. There is no unique way of 
choosing this parameterization. The simplest method comes in the form of Delves coordinates. Construction of these 
hyperangular coordinates is outlined in Appendix [A] and is described in detail in a number of references (see Refs. 
[911 197j for example) . This construction method also allows for a physically meaningful grouping of the cartesian 
coordinates. For example a hyperangular coordinate system that treats the dimer-atom-atom system as a separate 
three-body subsystem can be created. This type of physically meaningful coordinate system plays a crucial role in 
the construction of the variational basis set that follows. 

After adoption of the Jacobi vectors, the center of mass of the four-body system is removed, which leaves a 
9-dimensional partial differential equation to solve. By applying hyperspherical coordinates, this becomes an 8- 
dimensional hyperangular PDE that must be solved at each hyperradius, a daunting task. A further simplification 
is achieved by initially considering only zero total angular momentum states of the system. This implies that there 
is no dependence on the three Euler angles in the final wavefunction, and in a body-fixed coordinate system these 
three degrees of freedom can be removed. The body-fixed coordinates adopted here are called democratic coordinates, 
adequately described in several references (see Refs. [551 l9"51 199| ). The parameterization of Aquilanti and Cavalli is 
convenient for our purposes (for more detail see their work in Ref. [55] ). 

At the heart of democratic coordinates is a rotation from a space-fixed frame to a body-fixed frame: 



g = D T (a, 0, 7) g bf 



(26) 



where g is the matrix of "lab frame" Jacobi vectors defined in Eq. C12 g^f is the matrix of body-fixed Jacobi 



coordinates, and D (a, (3, 7) is an Euler rotation matrix defined in the standard way as 



D 



cos a — sin a 
sin a cos a 
1 



cos j3 sin 

1 
— sin j3 cos j3 



cos 7 — sin 7 
sin 7 cos 7 
1 



(27) 



This parameterization is described in detail in Appendix |A"[ After removing the Euler angles, the body fixed Jacobi 
vectors are then described by a set of 5 angles {0i, 02, $1, $2, ^3} and the hyperradius R. The angles ©i and ©2 
parameterize the overall x, y, and z spatial extent of the four-body system in the body-fixed frame, while the angles 
$1, $2, and $3 describe the internal configuration of the four particles. 

The description of coalescence points in democratic coordinates is especially important. These are the points at 
which interactions occur and also where nodes must be enforced for symmetry. Figure [2] shows these points for 
©i = it/2, which enforces planar configurations, and for several values of ©2. The body fixed coordinates in question 
are H-type Jacobi coordintes that connect identical fermions, so symmetry is is easily described. The $3 axis in Fig. 
[2] is shown from ir/2 to it and then to tt/2 to emphasize symmetry. The red surfaces surround Pauli exclusion nodes 
while the blue surfaces surround interaction points. It is clear that using a symmetry-based coordinate system leaves 
a simple description of the Pauli exclusion nodes. 



A. Unsymmetrized basis functions 



With the Jacobi vectors and democratic coordinates in hand, the 12-dimensional four-body problem is reduced to 
a 6-dimcnsional problem for total orbital angular momentum J = 0. After the hyperradius is treated adiabatically, 
the remaining 5-dimensional hyperangular partial differential equation Eq. (pi) must be solved at each R to obtain the 
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FIG. 2: Surfaces surrounding the coalescence points in the body-fixed democratic coordinates are shown for 9i = n/2 and 
82 = ^ (a), — (b), and — (c) respectively. Blue surfaces surround interaction coalescence points while red surfaces surround 
Pauli exclusion nodes. 

adiabatic channel functions and potentials used in the adiabatic hyperspherical representation. In Eq. V(R;H,) 
is chosen as a sum of short-range pairwise interactions, which to an excellent approximation affects only the s-wave 
for each pair: V(R;£l) = J2i j V ( r ij)> where the sum runs over all possible pairs of distinguishable fermions. This 
section only considers a potential whose zero energy s-wave scattering length a is positive and large compared with 
the range ro of the interaction. Further, unless otherwise stated, we assume that the potential can support only a 
single weakly-bound dimer. 

The strategy used here is not unknown [100] . It involves using a variational basis that diagonalizes the adiabatic 
Hamiltonian in two limits asymptotically (i? 3> a) and at small distances (R -C ro). It is thought that linear combi- 
nations of these basis elements will provide a variationally accurate description of the wavefunction at intermediate 
R- values. 

Next we describe the unsymmetrized basis functions that exactly diagonalize Eq. ^ in the small- R and large- R 
regimes. At large R, three scattering thresholds arise: a threshold energy corresponding to weakly-bound dimers at 
twice the dimer binding energy, another threshold consisting of a single weakly bound dimer and two free particles, 
and finally a threshold associated with four free particles. In general, it would be necessary to consider another set 
of thresholds associated with trimer states plus a free atom (for instance, a set of Efimov states for bosons) . But for 
equal mass fermions, such considerations are irrelevant since no weakly bound trimers occur in the a> ro regime. At 
small R, the physics is dominated by the kinetic energy, and the eigenstates of the adiabatic Hamiltonian are simply 
the 4-body hyperspherical harmonics which also describe four free particles at large R. For a detailed description of 
hyperspherical harmonics, see Appendix [B] Identification of these threshold regimes gives a simple interpretation of 
the corresponding channel functions and provides a starting point for the construction of our variational basis. 



1. Dimer- Atom- Atom Three-Body Basis Functions (2 + 1 + 1) 

One fragmentation possibility that must be incorporated into the asymptotic behavior of the four-fermion system is 
that of an s-wave dimer with two free particles. The dimer wave function 4>d is best incorporated using a hyperangular 
parameterization that treats the dimer-atom-atom system with a set of three-body hyperangles, described by 

*Wsb {Rj n) = ^ (ri2 ) y W3B ( n i3) ( (28) 

where Y\ 3Bfl3B is a three-body hyperspherical harmonic defined in Eq. ( |A13[ ), A 3 ^ is the three-body hyperangular 
momentum, and [133 indexes the degenerate states for each value of A3B. The dimer wave function tpd is chosen as 
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the bound state solution to the two-body Schrodinger equation: 



d 2 



2/i2fc Or 2 



V(r) 



rcj>d (r) = -E b r<pd (r) 



(29) 



Here the superscript 12 in fig^ indicates that the third particle in the three-body subsystem is a dimer of particles 1 
and 2. Further, for notational simplicity, [133 has been used to denote the set of quantum numbers, {l 2 , h, m 2 , 1713}, 
which enumerate the degenerate states for each X^b- 

So far the basis function defined by Eq. (28) is easily written in Delves coordinates. However, in order to ensure 



that L is a good quantum number, one must couple the angular momenta corresponding to the interaction Jacobi 
coordinates il (defined in Eq. CIO) to total angular momentum L = 0. The angular momentum of the (s-wave) dimer 
is by definition zero and all that remains is to restrict the angular momentum of the three-body sub-system to zero. 
This can be achieved by recognizing that the angular momentum associated with the individual Jacobi vectors are 
good quantum numbers for the hyperspherical harmonics defined by Eq. (A13), meaning that we can proceed using 
normal angular momentum coupling, i.e. 



*2+l+l 



(R,Q) = cb d (r 12 ) J2 E (hm 2 l 3 m 3 \00)Y X3Bf , 3B (nl%) 



(30) 



-l 2 m 3 = — 1 3 



where (limil-im^LM) is a Clebsch-Gordan coefficient, and l 2 (I3) is the angula r mo mentum quantum number as- 
sociated with p % 2 (P3 1 ) from the interaction Jacobi coordinates defined in Eqs. (CIO). Now with the total angular 
momentum set to L = 0, there must be no Euler angle dependence in the total wavefunction. The Delves coordinates 
can then be defined for this system in the body fixed frame using Eq. ( A9 ) . The Delves hyperangles are accordingly 



rewritten in terms of the democratic coordinates without including the Euler angle dependence. 



2. Four-Body Basis Functions + 1 + 1 + 

Another important asymptotic threshold that must be considered is that of four free particles. Using Delves 
coordinates, the free-particle eigenstates are four-body hyperspherical harmonics (see Appendix |b|) : 

(fi) =Nf? lmXl m Ntl mln sin A <- (<*,„,,«) cob 1 " (a, ro ,„) P^T-tyl (™s2a lm , n ) 

x N h j2 sin 1 ' (ai, m )cos /m {ai >m ) P lxt^-h-l m )/2 ( cos2a /,m) 
x Y limi (uji) Yi m „ lm (w m ) Yi nrHn (oj n ) , 

where fi has again been used to denote the set of quantum numbers {A12, li, h, ^3, tni, m<i, 7713} that enumerate the 
degenerate states for each A. Here ij is the spatial angular momentum quantum number associated with the Jacobi 
vector pi with z-projection m,, and A; >m is the sub-hyperangular momentum quantum number associated with the 



sub-hyperangular tree in Fig. 32 (For example, A12 = l\ + I2 + 2ni2 where 1212 is a non-negative integer.) The 
hyperangles {a; m; „, a^ m } are defined here using Delves coordinates as described in Appendix [A| and ui n refers to the 
spherical polar angles associated with the Jacobi vetor p n . 

The choice of quantum numbers described above does not give the total orbital angular momentum of the four 
particle system as a good quantum number. To accomplish this, the three angular momenta of the Jacobi vectors 
must be coupled to a resultant total , in this case to L = 0. This gives a variational basis function of the form 

*£?lt l ?{ty= £ E E E <W^ 3 m 3 |00) (31) 

Afi2 = — L12 TO3=— (3 JTl2 = — 12 mi = — h 

x (hmil 2 m 2 \L 12 M 12 ) (fi) . 

Now that the total angular momentum is set to L — the same procedure used for the vE^+i+i 2 basis functions can be 
employed. However, this time the hyperangular parameterization is defined using the symmetry Jacobi coordinates 



in Eqs. (C3). Since there is no dependence on the Euler angles, the Jacobi coordinates can then be defined in the 



body-fixed frame. 
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3. Dimer-Dimer Basis Functions (2 + 2) 

The asymptotic behavior of the two-component four-fermion system must include a description of two s-wave dimers 
separated by a large distance. To incorporate this behavior the variational basis must include a basis function of the 
form, 

* 2+2 (i?,ftH<Mr 12 )<Mr 34 ), (32) 

where the subscript 2+2 indicates the dimer-dimer nature of this function, and the dimer wavefunction, (f>d, is given 
by the two-body Schrodinger equation. Here /i 2 6 is the reduced mass of the two distinguishable fermions, and 



Eb « H /2/i26<2 is the binding energy of the weakly bound dimer. At first glance the right-hand side of Eq. (32) 



depends only implicitly on the hyperradius and hyperangles. To make this dependence explicit, Eqs. ( A22) and ( A27) 
are employed to extract r\ 2 (R, f2) and r 3 4 (R,fl). It can also be noted that the basis function, Eq. 32 does not 
respect the symmetry of the identical fermions, i.e. Pi 3 & 2 + 2 ^ —$2+2- The antisymmetrization of the variational 
basis is discussed in the next section. 

B. Symmetrizing the Variational Basis 

The definitions of the basis functions developed in the previous subsection do not include the fermionic symmetry 
of the four particle system in question. Until this point, we have only been concerned with Jacobi coordinate systems 
in which the particle exchange symmetry is well described and with a single set of Jacobi vectors that describe some 
of the interactions. In order to impose the S 2 <£> S 2 symmetry of two sets of two identical fermions, we now incorporate 
the extra H-type Jacobi coordinates described in Appendix [C] As a first step we define the projection operator, 

P=\(i-Pi 3 )(i-Pu), (33) 

where / is the identity operator, and Pij is the operator that permutes the coordinates of particles i and j. This operator 
will project any wavefunction onto the Hilbert space of wavefunctions that are antisymmetric under exchange of 
identical fermions. Since we are treating the fermionic species as distinguishable, permutations of members of different 
species are ignored. Applying this projection operator to the dimer-dimer basis function yields an unnormalized basis 
function, 

<T m) (R, fi) = P^2+2 (R, SI) = i {<j> d (r 13 ) <p d (r 34 ) - <i> d (r u ) 4> d (r 23 )) , (34) 



where the inter-particle distances r\± and r 2 3 are given in Eqs. (A24| and (A25) 



Imposition of the antisymmetry constraints on the dimer-atom-atom basis functions in Eq. ( 30 1 yields 



■y fa fa 

=^d(r ia ) E E (l2m 2 l 3 m 3 \00)Y X3Bfi3B (Ql 2 B ) (35) 

m 2 — — fa rns——l 3 
-i fa fa 

"4<Mr 23 ) E E (hm 2 l 3 m 3 \00)Y X3Btl3B (nf B ) 

1 h h 

- -^<l>d(ni) ^ Yl ( l 2m 2 l 3 m 3 \00) Y X3Btl3B 

m 2 =-h n»3=-l3 
-, h h 

+ J<l>d(ru) E E (l2m 2 l 3 m 3 \00)Y XsBIMsB (nl 



L 3B) ' 



m.2=— h rri3 = — h 



where Q 3B is the set of three-body hyperangles associated with particles i and j in a dimer and the remaining two 



particles free. The democratic parameterizations for the inter-particle distances from Eqs. (A22|-(A27| can be used 
in the dimer wavefunction directly. Through the use of symmetry coordinates, the hyperangles of the four-body 
system can be divided into a dimer subsystem and a three-body subsystem where the third particle is the dimer 
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itself. Using the three-body hyperangles in the three-body harmonic in each term of Eq. ( 35 1 , combined with the 



kinematic rotations from Eqs. (C14| and (C15), the three-body harmonics are then fully described in the hyperangles 



defined using symmetry Jacobi coordinates. Since 

^sy^n)X 3b l 2 l 3 been constramed tQ 

zero total spatial angular 

momentum, L = 0, the body-fixed parameterization of the Jacobi vectors can be inserted directly without worrying 
about the Euler angles a, (3 and 7. 

The final set of basis functions that must be antisymmetrized with respect to identical fermion exchange are the 
hyperspherical harmonics representing four free particles. Permutation of the identical fermions is accomplished in 



the symmetry coordinates using Eqs. (C4)-(C9). Using these permutations gives 



A 3 *iti+i+i («) = *i+i+i+i («) . 

where antisymmetry of the four free particle basis functions is enforced simply by choosing l\ and I2 to be odd. 

Another symmetry in this system is that of inversion (parity), in which all Jacobi coordinates are sent to their 
negatives, 

Pi -* 'Pi > 

where a = s,il,i2 and j = 1,2,3. Following the definitions of the Jacobi coordinates, positive inversion symme- 
try in the 1 + 1 + 1 + 1 basis functions, ^i+i+i+i (Q), is imposed by choosing A to be even. The 2 + 1 + 1 basis 

functions, ^jh^i+T^ 36 ' 2 ' 3 C^>^)> must already have positive inversion symmetry since <f>d (r) is an s-wave dimer wave- 
function and I2 = h f° r zero total spatial angular momentum, L = 0. The dimer-dimer basis function, ty^+T" 1 ^ (-^> 
is already symmetric under inversion and does not need further restrictions placed on it. 

The final symmetry to be imposed is not quite as obvious as the symmetries discussed so far. By performing a 
"spin-flip" operation in which the distinguishable species of fermions are exchanged, i.e. P12-P34, the Hamiltonian in 
Eq. [3] (with N = 4) remains unchanged. This operation is identical to inverting the two dimers in the dimer-dimer 
basis function. One can see that ty^+T™ 1 ^ 1S unchanged under this operation. We will limit ourselves to dimer-dimer 
collisions in this section and will only be concerned with basis functions that have this symmetry. This symmetry is 
imposed on both ^™™) A ^^ an d Vfigft* 1 ? by demanding l 3 to be even. 

Recalling that A = (li + V 2 + 2ni) + I3 + 2ro 2 where ri\ and n 2 are both non- negative integers, the combination of 
these symmetries implies that the minimum A for ^x+x+x+x mus t be X m in = 2. This argument plays a pivotal role 
in determining the overall threshold scaling law for four-body recombination, as is discussed in Section |IIB| 



IV. CORRELATED GAUSSIAN AND CORRELATED GAUSSIAN HYPERSPHERICAL METHOD 



A. Correlated Gaussian method 



In this Section, we discuss alternative numerical techniques to study the four-body problem. First, we present 
a powerful technique to describe few-body trapped systems where the solutions are expanded in correlated Gaus- 
sian (CG) basis set. Additional details regarding the CG basis set, including the evaluation of matrix elements, 
symmetrization, and basis set selection are discussed in Appendix [D] We then present an innovative method which 
combines the adiabatic hyperspherical representation with the CG basis set and Stochastic Variational method (SVM). 
For additional information on the methods described in this section, see ( |101[ 1102"] ). 



1. General procedure 



Different types of Gaussian basis functions have long been used in many different areas of physics. In particular, 
the usage of Gaussian basis functions is one of the key elements of the success of ab initio calculations in quantum 
chemistry. The idea of using an explicitly correlated Gaussian to solve quantum chemistry problems was introduced 
in 1960 by Boys |103j and Singer |104j . The combination of a Gaussian basis and the stochastical variational method 
SVM was first introduced by Kukulin and Krasnopol'sky [105] in nuclear physics and was extensively used by Suzuki 
and Varga [106H109] . These methods were also used to treat ultracold many-body Bose systems by Sorensen, Fedorov 
and Jensen |110) . A detailed discussion of both the SVM and CG methods can be found in a thesis by Sorensen |lllj 
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and, in particular, in Suzuki and Varga's book [112] . In the following, we present the CG method and its application 
to few-body trapped systems. 

Consider a set of coordinate vectors that describe the system {xi, ...,x N }. In this method, the eigenstates are 
expanded in a set of basis functions, 

*(cci, • • • ,x N ) = ®a{xi, ■ ■ ■ ,x N ) = ^2 C A (xi, ■ ■ ■ ,x N \A) . (36) 

A A 

Here A is a matrix with a set of parameters that characterize the basis function. In the second equality we have 
introduced a convenient ket notation. Solving the time-independent Schrodinger equation in this basis set reduces 
the problem to a diagonalization of the Hamiltonian matrix: 

UC, = EiOCi (37) 

Here, Ei are the energies of the eigenstates, Ci is a vector form with the coefficients Ca and H and O are matrices 
whose elements are Hba = (B\"H\A) and Oba = (B\A). For a 3D system, the evaluation of these matrix elements 
involves 3iV-dimensional integrations which are in general very expensive to compute. Therefore, the effectiveness of 
the basis set expansion method relies mainly on the appropriate selection of the basis functions. As we will see, the 
CG basis functions permit a fast evaluation of overlap and Hamiltonian matrix elements; they are flexible enough to 
correctly describe physical states. 

To reduce the dimensionality of the problem we take advantage of its symmetry properties. Since the interactions 
considered are spherically symmetric, the total angular momentum, J, is a good quantum number, and here we restrict 



ourselves to J = 0. Observe that if the basis functions only depend on the interparticle distances, then Eq. (36) only 
describes states with zero angular momentum and positive parity (J n = + ). Furthermore, in the problems we 
consider, the center-of-mass motion decouples from the system. Thus the CG basis functions take the form 




■■■ ,x N ) = i> Q {Rcu)S < exp - 2^ aijrfj/2 > , (38) 



where S is a symmetrization operator and nj is the interparticle distance between particles i and j. Here, ipo is 
the ground state of the center-of-mass motion. For trapped systems, Vo takes the form, iPo(Rcm) = e -fl cAf/ 2o >. . 
Because of its simple Gaussian form, ip can be absorbed into the exponential factor. Thus, in a more general way, 
the basis function can be written in terms of a matrix A that characterizes them, 

$ A (x 1 ,x 2l ...,x N ) = S \exp(-\x T .A.x) )■ = S i exp(-^ V AyXj.Xj) } , (39) 



2 



j exp(-^ £ Ai 



where x = {xi,x 2 , ...,xn} and A is a symmetric matrix. The matrix elements Ay are determined by the atj (see 



Appendix |D 3[ ). Because of the simplicity of the basis functions, Eq. (38), the matrix elements of the Hamiltonian can 
be calculated analytically. 

Analytical evaluation of the matrix elements is enabled by selecting the set of coordinates that simplifies the 



integrals. For basis functions of the form of Eq. (39), the matrix elements are characterized by a matrix M in the 
exponential. Hence the matrix element integrand can be greatly simplified if we write it in terms of the coordinate 
eigenvectors that diagonalize that matrix M. This change of coordinates permits, in many cases, the analytical 
evaluation of the matrix elements. The matrix elements are explicitly evaluated in Appendices |D 1| and |D 2] 

Two properties of the CG method are worth mentioning. First, the CG basis set is numerically linearly-dependent 
and over-complete, so a systematic increase in the number of basis functions will in principle converge to the exact 
eigenvalues |lllj . Secondly, the basis functions (5>a are square-integrable only if the matrix A is positive definite. 
We can further restrict the basis function by introducing real widths dij such that oiij = which ensures that 

A is positive definite. Furthermore, these widths are proportional to the mean interparticle distances in each basis 
function. Thus, it is easy to select them after considering the physical length scales relevant to the problem. Even 
though we have restricted the Hilbert space with this transformation, we have numerical evidence that that the results 
converge to the exact eigenvalues. 

The linear dependence in the basis set causes problems in the numerical diagonalization of the Hamiltonian matrix 



Eq. (37). Different ways to reduce or eliminate such problems are explained in the Appendix D 5 



Finally, we stress the importance of selecting an appropriate interaction potential. For the problems considered in 
this review, the interactions are expected to be characterized primarily by the scattering length, i.e., to be independent 
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of the shape of the potential. We capitalize on that flexibility by choosing a model potential that permits rapid 
evaluation of the matrix elements. A Gaussian form, 



V {r) = -dexp(-^j , (40) 

is particularly suitable for this basis set choice. If the range ro is much smaller than the scattering length, then the 
interactions are effectively characterized only by the scattering length. The scattering length is tuned by changing the 
strength of the interaction potential, d, while the range, ro, of the interaction potential remains unchanged. This is 
particularly convenient in this method since it implies that we only need to evaluate the matrix elements once and we 
can use them to solve the Schrddinger equation at any given potential strength ( or scattering length) . Of course, this 
procedure will give accurate results only if the basis set is complete enough to describe the different configurations 
that appear at different scattering lengths. 

In general, a simple version of this method includes four basic steps: generation of the basis set, evaluation of the 
matrix elements, elimination of the linear dependence, and evaluation of the spectrum. The stochastical variational 



where the basis functions are selected randomly. 



method (SVM), briefly discussed in Appendix D 6 combines the first three of these steps in an optimization procedure 



B. Correlated Gaussian Hyperspherical method 

Several techniques have been developed to solve few-body systems in the last few decades [38l I112H115] . Among 
these methods, the Correlated Gaussian (CG) technique presented in the previous Section has proven to be capable of 
describing trapped few-body systems with short-range interactions. Because of the simplicity of the matrix element 
calculation, the CG method provides an accurate description of the ground and excited states up to N = 6 parti- 
cles |116) . However, CG can only describe bound states. For this reason, it is numerically convenient to treat trapped 
systems where all the states are quantized. The CG cannot (without substantial modifications) describe states above 
the continuum nor the rich behavior of atomic collisions such as dissociation and recombination. 

The hyperspherical representation, on the other hand, provides an appropriate framework to treat the continuum. 
In the adiabatic hyperspherical representation (see Sec.|ll]), the Hamiltonian is solved as a function of the hyperradius 
R, reducing the many-body Schrodinger equation to a single variable form with a set of coupled effective potentials. 
The asymptotic behavior of the potentials and the channels describe different dissociation or fragmentation pathways, 
providing a suitable framework for analyzing collision physics. However, the standard hyperspherical methods expand 
the channel functions in B-splines or finite element basis functions [551 U17H119] , and the calculations become very 
computationally demanding for N > 3 systems. 

Ideally, we would like to combine the fast matrix element evaluation of the CG basis set with the capability of the 
hyperspherical framework to treat the continuum. Here, we explore how the CG basis set can be used within the 
adiabatic hyperspherical representation. We call the use of CG basis function to expand the channel functions in the 
hyperspherical framework the CG hyperspherical method (CGHS). 

In the hyperspherical framework, matrix elements of the Hamiltonian must be evaluated at fixed R. To proceed, 
consider first how the matrix element evaluation is carried out in the standard CG approach 

In the CG method, we select, for each matrix element evaluation, a set of coordinate vectors that simplifies the 
integration, i.e., the set of coordinate vectors that diagonalize the basis matrix M which characterizes the matrix 



element (see Appendix D 3 1 . The flexibility to choose the best set of coordinate vectors for each matrix element 
evaluation is key to the success of the CG method. 

The optimal set of coordinate vectors are formally selected by making an orthogonal transformation from an initial 
set of vectors x — {x%, x^} to a final set of vectors y = {yi, yjv}: xT = y, where T is the N x N orthogonal 
transformation matrix. The hyperspherical method is particularly suitable for such orthogonal transformations be- 
cause the hyperradius R is an invariant under them. Consider the hyperradius defined in terms of a set of mass-scaled 
Jacobi vectors [2U [321 [Ml 1120) . x = {xi, xat}, 

r 2 = T, x I ( 41 ) 

i 

If we apply an orthogonal transformation to a new set of vectors then 

r 2 = J2 x i = y TTjl y = E w? ( 42 ) 
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where we have used that T T T = I, and J is the identity. Therefore, in the hyperspherical framework we can 
also select the most convenient set of coordinate vectors for each matrix element evaluation. This is the key to 
reducing the dimensionality of the matrix clement integrals. One can view the flexibility afforded by such orthogonal 
transformations of the Jacobi vectors instead in terms of the hyperangles f2 that best simplify the evaluation of matrix 
elements. 

As an example of how the dimensionality of matrix-elements is reduced, consider a three dimensional iV-particle 
system in the center of mass frame and with zero orbital angular momentum ( J = 0). We will show that this technique 
reduces a (3N — 7)-dimensional numerical integral |121) to a sum over (N — 3) -dimensional numerical integrals (see 
Subsec. IV B 1). Hence, for N = 3 the matrix elements can be evaluated analytically, and the N = 4 matrix elements 
require a sum of one-dimensional numerical integrations. 

The next three subsections discuss the implementation of the CGHS. Many of the techniques used in the standard 
CG method can be directly used in the CGHS approach. For example, the selection and symmetrization of the basis 
function can be directly applied in the CGHS method. Also, the SVM method can be used to optimize the basis 
set at different values of the hyperradius R. Subsection |1V B 1| describes how the hyperangular Schrddinger equation 



(Eq. 43 ) can be solved using a CG basis set expansion and shows, as an example, how the unsymmetrized matrix 



elements can be calculated analytically for a four particle system. Finally, subsection IV B 2 discusses the general 
implementation of this method. 



1. Expansion of the channel functions in a CG basis set and calculation of matrix elements 

In the hyperspherical method (see Sec. [il]), channel functions are eigenfunctions of the adiabatic Hamiltonian 
H A (R;n), 

H A (R; n)<f>„(R; n) = U V {R)$ V {R; O). (43) 

The eigenvalues of this equation are the hyperspherical potential curves U V {R). The adiabatic Hamiltonian has the 
form, 

Here, d = 3Nj where At is the number of Jacobi vectors. 



A standard way to solve Eq. (431 is to expand the channel functions in a basis, 

\%(R;Q)) = J2^(R)\Bi(R;n)). (45) 



Here /x labels the channel functions and \Bj(R; fl)) are the CG basis functions [Eq. ( |39| )] written in hyperspherical 
coordinates. With this expansion, Eq. ( |43[ ) reduces to the generalized eigenvalue equation 

HAiRK = U^R)0(R) C)1 . (46) 

The vectors — {c^, cj^}, where D is the dimension of the basis set. Ha and O are the Hamiltonian and overlap 
matrices whose matrix elements are given by 

H A (Rh = ({BilHjL&tylBj)) , (47) 
0{R)n = ((BilBi)) . (48) 



Efficient evaluation of the matrix elements, e.g. Eqs. (47) and (48), is essential for the optimization of the ba- 



sis functions and the overall feasibility of the four-body calculations. Here, we demonstrate how to speed up the 
calculation by reducing the dimensionality of the numerical integrations involved in the matrix element evaluation. 

Consider a four-body system described by three Jacobi vectors, x = {xi, a; 2 , CC3}, once the center-of-mass motion 
is decoupled. The overlap matrix elements between two unsymmetrized basis functions $>a and <&b (characterized by 
matrices A and B in the respective exponents) is significantly simplified if we change variables to the set of coordinates 
that diagonalize A + B. We call /3±, ft% and fa the eigenvalues and y = {2/1,1/2, 2/3} are the eigenvectors of A + B. In 
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this new coordinate basis set the overlap integrand takes the form 



$ a (x 1 ,X2,x 3 )$b(xi,x 2 ,x 3 ,) = exp 



favl + favl + fay 2 



(49) 



In this set of eigencoordinates, the integration over the polar angles of t/i, vectors is easily carried out. To fix 
the hyperradius, we express the magnitude of the yi vectors in spherical coordinates, i.e. y\ — R sin 9 cos tf>, y 2 = 
i?sin6*sin(0)and y$ — Rcos9. In these coordinates the overlap matrix elements reads 



(B\A) „ = (4tt) 3 



exp 



R 2 {fa sin 2 9 cos 2 </> + fa sin 2 sin 2 + fa cos 2 9) 



>cos 2 ^>sin 2 0d6>#. (50) 



The integration over one of the angles can be carried out analytically. Introducing a variable dummy y, the overlap 
matrix clement takes the form 



(B\A) 



(47r) 3 7T 



R 2R 2 (fa-fa) 



exp 



R< 



[(fa+fa)(l-y 2 ) + 2fay 2 



h 



R- 



h-fa)(l-y 2 ) 



v 2 )dy, (51) 



where 1\ is the modified Bessel function of the first kind. 

To simplify the interaction matrix element evaluation, it is advantageous to use a Gaussian model potential as 
was used in the CG method. In this case, the interaction term can be evaluated in the same way as the overlap 

r 2 . 

term since the interaction is also a Gaussian. Each pairwise interaction can be written as Vij = Voexp(— ^p) = 

V exp(-x T .M^\x/(2d 2 )) (see Subsec." 
the structure 



D3 



for the definition of M^'). Therefore the interaction matrix element has 



(B\V ij \A)=V Jdtoexpi- 



x T .(A + B + M^/d 2 ).x 



(52) 



This integration can be performed following the same steps used for the overlap matrix element. Equation (|51|) can 



be used directly if we multiply it by Vb, and fa, fa and fa are replaced by the eigenvalues of A + B + M ^/djfNote 
that for each pairwise interaction, the matrix M^' changes and requires a new evaluation of the eigenvalues. 

The third term we need to evaluate is the hyperangular kinetic term at fixed R. This kinetic term is proportional 
to the grand angular momentum operator A, defined for the N — 3 case as 



A 2 fi 2 
VR 2 



Y> h 2 v 2 f^J_d_ R5 _d_ 

^ 2(i + 2(iR 5 dR OR' 



(53) 



The expression can be formally written as 



Tn =Tt — Tr, 



(54) 



where 



A 2 ^ 2 
2^lR 2 ' 



T t = - 



h 2 v 2 J ^ h 2 i d d5 d 

and Tr = p— i? 5 — . 

2/i ' H 2fiR 5 dR dR 



(55) 



In typical calculations, 7o is evaluated by directly applying the corresponding derivatives in the hyperangles O. 
However, in this case, it is convenient to evaluate 7r and Tr separately, since it is easier to differentiate over the 
Jacobi vectors and the hyperradius. These two matrix elements are not separately symmetric, but the angular kinetic 
energy matrix, i.e., the total kinetic energy minus the hyperradial kinetic energy, is symmetric. To obtain an explicitly 
symmetric operator, we symmetrize the operation (B\Tn\A) |^ = {{B\Tt ~ Tr|^4) \r + (A\Tt ~ Tr\B) \r)/2 and obtain 



(B\Tn\A) 



(4tt) 3 
i? 8 



exp 



favl + favi + favl 



T (y i , V2 , y3 )y i y 2 V3 d v 1 d V2dy, 



(56) 
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where 



3 r 



T(s/i,2/ 2) 2/3) = \ jlZ 



3ft + # - 2(A.B) U + 



R 2 



^ Ay? 



(y.A.y)(y.B.y) 



^ R 2 

\i=l 



B 2 



(57) 



It is easy to show that (A.-B)jj = Ylj=i a ijbij since A and B are symmetric matrices. Here the bar sign indicates the 
integration over the angular degrees of freedom of y±, y 2 , and y 3 . We then divide the total result by (47r) 3 . Making 
these integrations analytically we obtain 



(y.A.y) {y.B.y)=Y^ a u b H yf + ^ 



i=l 



i>3 



Hi bjj 



"<A, J ViVy 



(58) 



Rewriting the variables in spherical coordinates, we separate the hyperradial dependence in Eq. (56 1. As in Eq. (50 1, 



one of the angular integrations can be evaluated analytically and the final expression reduces to a one dimensional 
integral involving modified Bessel function of the first kind (see Ref. |101j for more details). 

The matrix elements involved in the B and Q couplings can be evaluated by following the above strategy, and it 
also reduces to a one dimensional numerical integration. The symmetrization of the matrix elements is handled just 
as in the standard CG method and is described in Appendix |D 1| 



2. General considerations 



Many of the procedures of the standard CG method can be easily extended to the CGHS. The selection, sym- 



Dl 


D3 


D4 


D5 



and D 6 1. However, the evaluation of the unsymmetrized matrix elements at fixed B is clearly different. Furthermore, 
the hyperangular Hamiltonian [Eq. 43 needs to solved at different hyperradii B. 

There are several properties that make the CGHS method particularly efficient. For the model potential used, the 
scattering length is tuned by varying the potential depths of the two-body interaction. Therefore, as in the CG case, 
the matrix elements need only be calculated once; then they can be used for a wide range of scattering lengths. Of 
course, the basis set should be sufficiently complete to describe the relevant potential curves at all desired scattering 
length values. 

The selection of the basis function generally depends on B. To avoid numerical problems, the mean hyperradius of 
each basis function (R) B should be of the same order of the hyperradius R in which the matrix elements are evaluated. 
We can ensure that (R) B ~ R by selecting some (or all) the weights d^ to be of the order of R. 

We consider two different optimization procedures. The first possible optimization procedure is the following: First, 
we select a few basis functions and optimized them to describe the lowest few hyperspherical harmonics. The widths 
of these basis functions are rescaled by R at each hyperradius so that they represent the hyperspherical harmonics 
equally well at different hyperradii. These basis functions are used at all R, while the remaining are optimized at 
each R. Starting from small R (of the order of the range of the potential), we optimize a set of basis functions. As R 
is increased, the basis set is increased and reoptimized. At every R step, only a fraction of the basis set is optimized, 
and those basis functions are selected randomly. After several i?-steps, the basis set is increased. 

Instead of optimizing the basis set at each R, one can alternatively try to create a complete basis set at large 
Rmax- ln this case, the basis functions should be complete enough to describe the lowest channel functions with 
interparticle distances varying from interaction range ro up to the hyperradius R max ■ Such a basis set can be rescaled 
to any R < R m ax and should efficiently describe the channel functions at that R. The rescaling procedure is simply 
dij I R = d™ ax / R max . This procedure avoids the optimization at each R. Furthermore, the kinetic, overlap, and 
couplings matrix elements at R are straightforwardly related to the ones at R ma x- The interaction potential is the 
only matrix element that needs to be recalculated at each R. This property can be understood by dimensional 
analysis. The kinetic, overlap and couplings matrix elements only have a single length scale R, so a rescaling of the 
widths is simply related to a rescaling of the matrix elements. In contrast, the interaction potential introduces a new 
length scale, so the matrix elements depend on both R and d , and the rescaling does not work. 

These two choices, the complete basis set or the small optimized basis set, can be appropriate in different circum- 
stances. If a large number of channels are needed, the complete basis method is often the best choice. But, if only a 
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small number of channels are needed, then the optimized basis set might be more efficient. 

The most convenient method we have found to optimize the basis functions in the four-boson and four-fermion 
problem is the following: First we select a hyperradius R m that is R m s» 300 do where the basis function will be 
initially optimized. The basis set is increased and optimized until the relevant potential curves are converged and, 
in that sense, the basis is complete. This basis is then rescaled, as proposed in the second optimization method, 
to all R < R m . For R > R m , it is too expensive to have a "complete" basis set. For that reason, we use the first 
optimization method to find a reliable description of the lowest potential curves. 

Note that for standard correlated Gaussian calculations, the matrices A and B need to be positive definite. This 
condition restricts the Hilbert space to exponentially decaying functions. In the hyperspherical treatment, this is not 
necessary since the matrix elements are always calculated at fixed R, even for exponentially growing functions. This 
gives more flexibility in the choice of optimal basis functions. 



V. APPLICATION TO THE FOUR-FERMION PROBLEM 



This section presents our results for the four-body fermionic problem using the methods discussed in Sect ions [TTT] and 



IV Our finite-energy calculations for elastic and inelastic processes are compared to established zero-energy results 
and are seen to exhibit significant qualitative and quantitative differences. Several properties of trapped four-fermion 
systems are also discussed, along with the connections between this few-body system and the many-body BEC-BCS 
crossover physics. 



A. Four-fermion potentials and the dimer-dimer wavefunction 



Calculation of the hyperradial potentials and channel functions using the variational basis method of Sec. Ill is 
conceptually simple. Matrix elements of the hyperangular part of the full Hamiltonian are required, 

h 2 A 2 . 

where the sum runs over all interacting pairs of distinguishable fcrmions. Sec. |III| considered the specific two-body 
interaction to be general, but required the two-body potential to support a weakly bound dimer state (and hence 
a positive scattering length much larger than the range of the interaction). At this point we adopt the so-called 
Poschl- Teller potential, 

V(r) = -g— , (59) 

cosh (r/ro) 

where ro is the range of the interaction. Unless otherwise stated Uq is tuned so that V (r) gives the appropriate scat- 
tering length with only a single shallow bound state. This potential is adopted because the bound state wavefunctions 
and binding energies are known analytically [122], but any two-body interaction could be used here, provided that 
one obtains the wavefunctions and energies numerically or analytically. 

Application of the variational basis results in a generalized eigenvalue problem, 

H (R) x v (R) = U v (R) S [R) x v (R) (60) 

where U v (R) is the v-th adiabatic hyperradial potential, and x v is the channel function expansion in the variational 
basis. The matrix elements of H are given by matrix elements of the adiabatic Hamiltonian at fixed hyperradius, 

Because the variational basis is not orthogonal, a real, symmetric overlap matrix, S, appears in this matrix equation. 
While the method employed here is conceptually simple, the actual calculation of the matrix elements is numerically 
demanding because the interaction valleys in the hyperangular potential surface, - V (fij), become localized into 
narrow cuts of the hyperangular space at large hyperradii. Further, examination of Fig. [2] shows that the locus of 
coalescence points where the interatomic potential is appreciable has a complicated structure in the five dimensional 



body-fixed hyperangular space. To accurately calculate the matrix elements in Eq. (60 1 numerically, a large number 
of integration points must be placed within the interaction valleys. 
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Despite all of these complications, the adiabatic potential can be found approximately. Figure [3] shows the full set 
of hyperradial potentials including the diagonal non-adiabatic correction (solid curves) calculated using 8 variational 
basis elements: one 2 + 2 element, four 2 + 1 + 1 elements, and three 1 + 1 + 1 + 1 elements. Also shown are the 
expectation values of the basis elements themselves, i.e. the diagonal of H(R) from Eq. (60) (dashed curves). All 
calculations shown here are performed for a = lOOro It is clear that the lowest potential curve converges very quickly 
with respect to the number of variational basis elements used. The lowest potentials converge well when only a few 
variational basis functions are included, while the higher potentials are somewhat suspect. According to the universal 
theory of zero-range interactions, the hyperspherical potential curves should only depend on a in the regime where o 
is the dominant length scale in the problem. Thus, in our finite range interaction calculations, the adiabatic potentials 
should become universal in the i? 3> tq regime for large scattering lengths, i.e. a>ro. In other words, the potentials 
should look the same when scaled by the scattering length and the binding energy, U V (R 3> r ) = (fi 2 /ma 2 )u u (R/a) 
where u v (x) is a universal function for the z/th effective potential. Comparison with the potential curves computed 
in the correlated Gaussian method shows excellent agreement in the lowest dimer-dimer potential, and reasonable 
agreement for the lowest few dimer-atom-atom potentials jlOlj . 




FIG. 3: The calculated hyperradial potentials (solid lines) for a = lOOro are shown as a function of R/a. Also shown are the 
expectation values of the fixed- R Hamiltonian for the individual variational basis functions (dashed curves). 

At large R, the lowest hyperradial adiabatic potential curve (see Fig. [3| approaches the bound-state energy of two 
dimers that are approximately separated by a distance R. It is natural to interpret processes for which flux enters 
and leaves this channel as "dimer-dimer" collisions. Examining this potential further, one can see that at hyperradii 
less than the scattering length, R < a, the adiabatic dimer-dimer potential becomes strongly repulsive. This can be 
visualized qualitatively as hard wall scattering, which would give a dimer-dimer scattering length comparable to the 
two-body scattering length ~ a. Higher potential curves approach the single dimer binding energy at large R, 
indicating that these potentials correspond to a dimer with two free particles in the large R limit. Note that the 



variational basis functions described in Section [HI] give the correct large R adiabatic energies by construction. As the 
scattering length becomes much larger than the range of the two-body potential, the effective four-fermion hyperradial 
potential becomes universal and independent of a. In the range of r <C R <C a: 

2+1 ' (61) 

where po = 2.55. This universal potential was extracted in Refs. |123| i!24j by examining the behavior of the ground 
state energy of four fermions in a trap in the unitarity limit. 

Figure |4] shows the coupling strengths, h 2 P 2 m / {2fi \(U m (R) — U n {R))]}, between the dimer-dimer potential and 
the lowest three dimer-atom-atom adiabatic potentials for a two-body scattering length of a — lOOro. In each case the 
coupling strength peaks strongly near the short range region, R ~ r , and near the scattering length, R ~ a, and then 
falls off quickly in the large R limit. This behavior indicates that recombination - from a state consisting of a dimer 
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and two free particles to the dimer-dimer state - occurs mainly at hyperradii of the order of a. Looking at Fig. [4] 
one might think that a recombination path which occurs at small R, R ~ ro, could also contribute, but the strong 
repulsion in the dimer-atom-atom potentials between R ~ r and R ~ a, shown in Fig. [3j suppresses this pathway. 

0.005 
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FIG. 4: The nonadiabatic coupling strength between the dimer-dimer potential and the lowest dimer-atom-atom potential is 
shown as a function of R (ro = 100 a.u. was choosen to be the van der Waals length of 4 0K). The blue dashed line shows the 
position of the coupling peak at R/a fn 3.5. 

Figure [5] shows an isosurface of the hyperangular probability density in the conhgurational angles {<fti, <j>2: 03 } after 
integrating out 0i and 02 at a fixed hyperradius of R = 0.41a. The <f>i axis has been modified here by shifting the 
region n/2 < (f>i < n to emphasize the symmetry of the system. Each cobra-like surface corresponds to a peak in the 
four-body probability density. By examining Fig. [2] one sees that the spine of each cobra corresponds to the locus of 
interaction coalescence points. For each choice of {4>i, §2-, </>3}j the maximum of the probability density in &i and 2 
is given in a planar geometry, Oi = tt/2. The coloring of each cobra indicates the value of O2 at which the maximum 
occurs. Darker colors indicate a more linear geometry, i.e. 02 is closer to 0. Figure [6] shows the same plot for the 2 + 2 
basis function only. A comparison of Figs. [5] and [6] indicates that the added variational basis elements are critical for 
describing the full dimer-dimer channel function for hyperradii less than the scattering length. 



B. Elastic dimer-dimer scattering 

With the hyperradial potentials and non-adiabatic couplings in hand, low energy dimer-dimer scattering properties 
can be examined. The zero-energy dimer-dimer scattering length in the limit of large two-body scattering length was 
first calculated by Petrov et. al [125J and found to be 

a dd (0) = 0.60 (2) a, (62) 

where the number in the parentheses indicates ±0.02, the 2% error stated in Ref. 125 . This result has been confirmed 
using several different theoretical approaches [751 11231 1124L 1126) . 

Using the adiabatic potentials shown in Fig. [3] and the resulting non-adiabatic couplings, the energy-dependent 
dimer-dimer scattering length defined by 

add (E co \) = : (63) 

kdd 

can be calculated. Here E co \ is the collision energy of the two dimers with respect to the dimer-dimer threshold, 
and Sdd is the s-wave dimer-dimer phase shift. When the collision energy becomes greater than the dimer binding 
energy, the two dimers collide with enough energy to potentially dissociate one of them. When this happens, the 
four fermion system can fragment in an excited hyperspherical channel causing a loss of flux from the dimer-dimer 
channel. This process is parameterized by the imaginary part of the dimer-dimer scattering length becomes non-zero 
when E co \ > Et,. 
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FIG. 5: An isosurface of the dimer-dimer probability density is shown. The surfaces are found by integrating the total 
probability over 9\ and #2 and plotting with respect to the remaining democratic angles (<f>i, cj>2, <j>3)- The peak probability 
always occurs in planar configurations, 9\ = 7r/2. The coloring (light to dark) indicates the value of 82 at the peak. 



Figures [7] and [8] respectively show the real and imaginary parts of the dimer-dimer scattering length calculated 
with different numbers of adiabatic channels plotted as functions of E co \ in units of the dimer binding energy. Also 
shown in Fig. [7] is the dimer-dimer scattering length calculated from the variational potential that results from using 
a single variational basis element. It is important to note that the single adiabatic channel calculation and the single 
basis function calculation are not the same. In the former, the single potential used is the lowest potential resulting 
from a calculation using multiple basis functions, while the latter is the result of using only the 2 + 2 variational basis 
function and is guaranteed to be less accurate. Not surprisingly, the scattering length at collision energies comparable 
to the binding energy depends strongly on the number of channels used. With just a single channel in use, there is 
no decay pathway available for the system. As more channels are included the system has more pathways into which 
it can fragment, which modifies the high energy behavior. 

What is more surprising is the low energy behavior seen in Fig. [7] For a single variational basis element, the 
dimer-dimer zero energy scattering length is found to be add = 0.72a, which is already within 20% of the result of 
Ref. |125) , add (0) — 0.6a. A single channel calculation using the dimer-dimer potential and channel function that 
results from using 5 basis elements improves considerably on this yielding add (0) = 0.64a, showing that inclusion 
of correlations characteristic of two free particles at hyperradii less than a gives a significant contribution to the 
physics of dimer-dimer scattering. It is somewhat unexpected that the single channel calculation is only 8% off of 
the predicted value. As the scattering energy approaches zero, the higher fragmentation channels become strongly 
closed but still apparently play a small role in the dimer-dimer scattering process. By including progressively more 
channels in the scattering calculation the zero-energy dimer-dimer scattering length can be extracted for large two 
body scattering length: 



a dd (0) = 0.605 (5) a. 



(64) 



This result is in agreement with the results of Ref. [123LI124] and the results of Section V F which found the zero-energy 
dimer-dimer scattering length to similar accuracy using different methods. 



C. Energy Dependent Dimer-Dimer Scattering 



By examining the low energy behavior of the energy dependent dimer-dimer scattering length, the effective range can 
be extracted. The two dimers "see" each other when their wavefunctions are overlapping, i.e. when the hyperradius 
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FIG. 6: The same as Fig. [5] but using only the 2 + 2 function. The dashed gray lines are purely for perspective. 
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FIG. 7: The real part of the energy dependent dimer-dimer scattering length is shown in units of the atom-atom scattering 
length a, plotted versus the collision energy in units of the dimer binding energy. The calculation is carried out with one, two, 
three, four, and five adiabatic channels (blue, black, red, green, and purple curves respectively) from the adiabatic potential 
curves and couplings computed using 8 basis functions. The red dashed line shows add = 0.6a, the prediction of Ref. [127] . 



is approximately equal to the scattering length, R ~ a. If one thinks of the effective range of an interaction as 
proportional to the size of the interaction region, then one would expect the effective range for dimer-dimer scattering 
to be proportional to the scattering length. By fitting the low energy scattering phase shift to the effective range 
expansion, 

t 1 — = k dd cot S d d = -7777 + \r dd k dd 2 + 0(k dd 4 ) (65) 

a dd (E) a dd (0) 2 

this intuitive behavior is born out, giving an effective range: 



r dd = 0.13(l)a, 



(66) 
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FIG. 8: The imaginary part of the energy dependent dimer-dimer scattering length is shown in dimensions of a plotted versus 
the collision energy in units of the binding energy. The calculation is carried out with one, two, three, four, and five adiabatic 
channels (blue, black, red, green, and purple curves respectively) from the adiabatic potential curves and couplings computed 
using 8 basis functions. 



where a is the two-body scattering length. Figure [9] shows both the real and imaginary parts of the energy dependent 
dimer-dimer scattering length as a function of collision energy in units of the binding energy compared to the effective 
range expansion, Eq. (65 1. This clearly shows that, while the low energy behavior of dimer-dimer scattering is 



well described by the effective range expansion, it is only accurate over a small range of collision energies. In fact, 
for collision energies larger than the binding energy, add (E co {) actually turns over and decreases as dimer breakup 
channels become open. Further, when the collision energy exceeds the dimer binding energy, E co \ — E^, the dimer- 
dimer scattering length becomes complex, with an imaginary part that parameterizes inelastic processes. These results 
indicate that both the real and imaginary dimer-dimer scattering lengths are universal functions of the collision energy. 
Specifically, they are insensitive to the short range nature of the two-body interaction, for scattering lengths much 
larger than the two-body interaction length scale, r . Because very few basis functions were used in these calculations, 
the results at higher energies, Ef, <C E co \ <C h 2 /mr 2 , are n °t weu converged, though their qualitative nature is expected 
to persist. Well above the dissociation threshold, the real part of add exhibits an oscillatory behavior. These oscillations 
are caused by interference between different scattering pathways. As more basis functions are included and the high 
energy results converge, the large number of available pathways generally washes out the oscillatory behavior and 
produces incoherence , but the decrease in the real part of the add at higher energies is expected to survive as the 
calculations become better converged. 

The dependence of add on a at finite collision energy is particularly interesting. In the large a limit, the dimer 
binding energy becomes Ej, w h 2 /ma 2 , so that as a increases, the binding energy decreases. At the critical value of 
the scattering length, 



\fmEcoi ' 

the collision energy coincides with the binding energy, and the dimer-atom-atom channel becomes open. As a result 
(see Fig. [9]), one expects the real part of add to turn over and remain finite for all values of a. This behavior is 
demonstrated in Fig. 10 which compares the real part of the dimer-dimer scattering length at several fixed collision 
energies with the zero-energy result, add (0) = 0.6a. The scattering length scale has been fixed by setting the range 
of the interaction to be approximately the Van der Waals length of 40 K, r sa 100 a.u. Another aspect of the finite 
collision energy behavior is that at large scattering length, the dimer-atom-atom channels become open, and dimer 
dissociation is allowed. Thus, near unitarity the Fermi gas might be viewed as a coherent mixture of atoms and weakly 
bound dimers. 



D. Dimer-dimer relaxation 



A significant loss process in an ultracold gas of bosonic dimers is that of dimer-dimer relaxation, in which two 
dimers collide and in the process at least one of the dimers relaxes to a deeply bound two-body state. The extra 
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FIG. 9: The real (red) and imaginary (green) parts of the energy dependent dimer-dimer scattering length are shown in units 
of a, plotted versus the collision energy in units of the binding energy. Also shown is the energy dependent scattering length 
predicted by the effective range expansion. Figure from Ref. |75j . 
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FIG. 10: The real part of the energy dependent dimer-dimer scattering length is shown as a function of the two-body scattering 
length in atomic units for several collision energies: E^/kb =250nK,100nK, 25nK, lOnK, 2.5nK, InK, lO^nK, and 10 _2 nK. 
Also shown is the zero energy prediction (black dashed curve). Figure from Ref. |75j . 



binding energy is released as kinetic energy which is sufficient to eject the remaining fragments from the trap. This 
process was studied by Petrov, Salomon and Shlyapnikov [1251 1127j , who assumed that the relaxation rate is controlled 
by the probability for three particles to be found in close proximity to one another. With this assumption and the 
further assumption that the fourth particle is far away and plays no role in the scattering process, they predict that 
the relaxation rate is suppressed at large two-body scattering lengths with a scaling law, V~i oc a -2 ' 55 . 

Here we introduce a new method for finding the dimer-dimer relaxation rate based directly on Fermi's golden rule. 
The key observation in this section is that the final allowed states appear as an infinite set of hyperspherical potentials 
corresponding to a deeply bound dimer with two free atoms. The transition rate to a single one of these potentials 
can be described by the Fermi-popularized golden rule, i.e. 



Tp <x | (* dd {R; Si) \V (R, SI) | * A (R, Sl))\ 2 



(67) 
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Here is the final outgoing state, ^dd is the dimer-dimer wavefunction, and V (R, ft) is the sum of the two-body 
interactions. This matrix element and the sum of probabilities over final states are evaluated in Appendix [E] The 
final result of this analysis is expressed as an integral over the hyperradius, 



PWKB (R)J'(R) 

Rk (R) 



p (R) dR 



(68) 



where Py/K b (R) is the WKB probability density of the dimer-dimer wavefunction at hyperradius R, k (R) is the 
WKB wavenumber: 
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(69) 



In Eq. (68) p(R) is the nearly constant density of final states, and T (R) is the probability for three particles to be 



near one another in the dimer-dimer wavefunction at hyperradius R: 

T(R) = (R;fl)\f(R,n)\^ dd (iJjfi)) 



(70) 



Here &dd is the hyperangular dimer-dimer channel function, and / (i?, Q) is a proximity function that is appreciable 
only when three particles are all approximately within the range of the two-body interaction. 



Equation ( 68 ) makes physical sense upon closer examination. It says that the rate at which a dimer relaxes to a 



deeper state is determined, with some extra factors, by the probability that three particles are close enough together 
so that two of them can fall into a deeply bound state and release the extra binding energy to the third particle. 
Figure 11 shows the integrand from Eq. (68 1 for several scattering lengths as a function of the hyperradius in units 



of the scattering length. This quantity can be interpreted as being proportional to the transition rate per unit 
hyperradius, i.e. the probability that the transition will occur between R and R + dR. The full transition rate is 
determined by the nature of the interaction at short range and is not predictable using this method. By examining the 
relaxation rate as a function of scattering length, however, a scaling law can be extracted at each fixed hyperradius. 




FIG. 11: The integrand from Eq. (681 is shown for a = 50ro (red),64ro (green), 80ro (blue), and lOOro (black) as a function of 
R/a. 



Figure 12 shows the relaxation rate per unit hyperradius for several fixed values of R/a as a function of the scattering 
length, a. The large a behavior in each case appears to follow a scaling law, but the scaling law changes with R/a. 



This behavior indicates that, contrary to the prediction of Ref. |127j . when the integral in Eq. (68) is evaluated, 



the relaxation rate will not be determined by a simple power law. By integrating over different hyperradial regions, 
contributions to the transition rate from different processes can be extracted. For instance, if the integral in Eq. (68) 



is performed only over small hyperradii, R < 5rg, the result is the transition rate due to processes in which all four 
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FIG. 13: The relaxation rate in a dimer-dimer collision is shown as a function of the atom-atom scattering length (see text). 
Figure from Ref. [75] 

Figure [13] shows the relaxation rate as a function of the scattering length in atomic units as a solid line. In this 
result the range of the interaction is set to the van der Waals length of 40 K, tq ss 100 a.u. Also shown in Fig. 13 



are the contributions to this relaxation rate due to four-body processes (dashed blue curve), and due to three-body 
processes (dotted green curve). Also shown is the expected scaling law for transitions that occur at small hyperradius, 
R = 5rQ. Because the hyperradius is small in this regime, the probability of three particles being in proximity is near 
unity, meaning that the transition probability per unit hyperradius is determined by the probability that the system 
can tunnel through the repulsive potential seen in Fig. 3' at R < a. The universal repulsive potential in this regime 

~ 2/i ' (71) 

Po = 2.55, (72) 



leads to a scaling law for transitions in the small R regime that behaves as 



(73) 



Figure 13 also shows the experimentally determined relaxation rates from Ref. 



Both the scaling law predicted in 
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Ref. |127j of a~ 2 - 55 and the prediction using Eq. (68) are consistent with the experimental data in the regime for 1000 
a.u.< a < 4000 a.u.. The experimental data for a > 3000 a.u. are in the regime where the average dimer separation 
is less than the dimer size, where the dimer-dimer scattering picture discussed here no longer applies. 

E. Trapped four-body system 

In this section we abandon the hyperspherical methods of the previous sections and examine the case of trapped 
four-fermion systems. The four-body system in confined geometries has recently become computationally accessible. 
In particular, the trapped two-component Fermi system has been intensely studied in the last few-years and has 
become a benchmark for different theories and universal predictions. One of the challenges present in the study 
of universal four-body physics in dilute gases with short range forces is the description of disparate length scales 
associated with the large interpaticle distances and the short-range two-body interactions. 

In this section, we analyze the spectrum, dynamics and universal properties of four-body solutions in confined 
geometries. The four-body system is described by a the model Hamiltonian 

" = 1 + r«^) + 1 (ss*J + ?^ 2 "?) + 1 1 nr,) (74) 

i— 1 v 7 i' — l v 7 i— 1 i' = l 

where unprimed indices label the fcrmionic species with mass mi, primed indices label the species with mass m.2, 
and Ti is the position vector of the ith fermion. The trapping frequency u> is assumed to be equal for both species. 



In order to facilitate a calculation with the CG method described in Section IV we take the interaction potential 



V to be a purely attractive Gaussian (see Eq. (40)) and tune the depth of V to give the desired (large) scattering 



length. The mass ratio k is defined by mi/mj, and throughout the analysis we assume mi > m-z. A trap length 
a^l — y^h/rriiUj is defined for each species as well as a trap length associated with the pair a^o — ■\Zh/2/jw, where 

fi = mi 7712/ (mi + JTI2). For equal mass systems, a^o = = a?„- 

This section reviews a series of predictions for the two-component four-fermion system. The spectrum and structural 
properties of the four-fermion system are analyzed throughout the BCS-BEC crossover, followed by an exploration of 
the system dynamics as the scattering length is tuned close to a Fano-Feshbach resonance. Finally, we review a series 
of numerical studies that confirm and quantify universal predictions. 

1. Spectrum in the BCS-BEC crossover 

To obtain the J = spectra for the four-fermion system in the BCS-BEC crossover problem, the CG method 
(Section |IV| is utilized to solve the time-independent Schrodinger equation for different values of a. Like most 
numerical methods, this method provides an adiabatic spectrum (in time), i.e., the energies of the spectrum are 
labeled according to their energy values as a changes. The four-body spectra present a series of crossings or narrow 
avoided crossings when the scattering length is tuned across the BCS-BEC crossover. For this reason, it is convenient 
to use a representation where these narrow avoided crossings are treated diabatically, and the spectrum smoothly 
evolves from the BCS to the BEC side. The diabatic representation is more relevant from the physical point of view 
since the diabatic states are usually associated with good or "approximately good" symmetries of the problem. 

To illustrate the diabatization procedure, consider the spectrum of the four-fermion system in the strongly interact- 



ing region shown in Fig. 14 A series of crossings and avoided crossings occurs when the adiabatic parameter A = 1/a 
is varied in the strongly interacting region. The avoided crossings can be roughly characterized by their width AA, 
the range where the two adiabatic energy curves are coupled. There are narrow crossings where A A <C l/a?io and 
there are wide crossings where AA > 1/aho- To obtain smooth energy values, we use variation of the diabatization 
procedure presented in Ref. |128j . 

The objective of the diabatization algorithm is to make the one-to-one connection between states and energies in 
consecutive points of the A grid that maximize the sum of the overlaps between connected states. The diabatization 
procedure starts from the BCS (a < 0) side of the resonance and connects the states (and their energies) between 
consecutive values of A for which their overlap is maximum. When two initial energies connect to the same final 
energy, a refinement of the diabatization procedure is applied. 

Diabatization is controlled by the spacing between consecutive values of A given by AA 9 . If the width of the avoided 
crossing is smaller than AA g , then that crossing is diabatized. But if the width of the avoided crossing is larger than 
AA g , then that crossing is not diabatized. Thus, AA g is selected so that narrow crossings are diabatized and wide 
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FIG. 14: Four-fermion energy spectrum as a function of A = 1/a in the unitarity region with J = 0. The thin solid black lines 
correspond to the adiabatic spectrum. The wide black line with circles is the diabatic ground state labeled $cm- The blue 
curve with circles is the diabatic first-excited state labeled $ dd2, the wide red curve with circles is the diabatic state ^daa, 
and the wide green curve with circles is the diabatic state ^4A- 



crossings remain adiabatic. For example, in Fig.[l4]we see how this procedure diabatizes the narrow crossings of 
however, wide crossings such as the one between \1/ daa and ^aa are still adiabatic in this representation. 

This structure of avoided crossings permits a global view of the manner in which states evolve from weakly interacting 
fermions at a < to all the different configurations of a Fermi gas at a > 0, i.e., molecular bosonic states, fermionic 
states, and molecular Bose-Fermi mixtures. Furthermore, it allow us to visualize concretely the alternative pathways 
of the time-dependent sweep experiments. 



The diabatic spectrum of the four-fermion system is presented in Fig. 15 The structure of avoided crossings is 
complicated because of two different thresholds exist, one corresponding to the dimer-atom-atom and one to the 
dimer-dimer state. We identify three different families of diabatic states in this spectrum. The dimer-dimer family, 
represented by the black and blue curves, describes the ground and excited dimer-dimer states. These states are 
separated by approximately 2Huj on the BEC side. The dimer-two-atom family, represented by the red energy curves, 
follows the dimer binding energy. In the BEC limit, the dimer-two-atom family reproduces the degeneracies of three 
distinguishable particles: a spin up atom, a spin down atom, and a dimer. The third family describes four-atom 
bound states, for which the atoms form no dimers, and whose energy remains positive in the crossover region. In the 
BEC limit as a — > 0, the four-atom family reproduces the spectrum of the nonintcracting four-body system. 

The evolution of the N — 4 spectra through the BCS-BEC crossover region can be understood qualitatively by 
considering the important quantum numbers for the description of the dimer. For each vibrational excitation of 2huj 
in the noninteracting limit, there is one state that diabatically becomes a dimer-dimer state. These states correspond 
qualitatively to states where the relative angular momentum of two spin- up-spin-down pairs is zero [L^ el — 0] , and the 
relative angular momentum between the pairs is also zero. The spin-up-spin-down pairs are in the lowest vibrational 
state. In the weakly interacting BCS limit, where the degeneracy of the vibrational states is broken, pair-pair states 
correspond to the lowest states. 

A direct and more concrete way to visualize the structure of the spectrum is to analyze the evolution of the adiabatic 
hyperspherical potential curves. Figure 16 presents the four-fermion adiabatic hyperspherical potential curves U V {R), 
obtained with the correlated Gaussian hyperspherical method (CGHS). Panel (a) presents the potential curves in the 
BCS regime, which are clearly grouped into families. Potential curves belonging to the same family are degenerate 
in the noninteracting limit. Thus, the weak interactions in the BCS regime break the degeneracies of the potential 
curves forming these families of potential curves. Panel (b) describes the system in the BEC regime. In this case, 
the description of the system is quite clear. The lowest potential curve is more than twice as deep as the rest of the 
curves and is associated with the dimer-dimer threshold. The family of dimer-dimer states live mainly in the lowest 
potential curve. The remainder of the displayed potential curves are associated with the dimer-two-atom threshold. 
The dimer-two-atom states are mainly described by this family of potential curves. A third family of potential curves, 
not shown in Fig. 16 (b), describes four-atom states. This family of potential curves has a different large- R asymptotic 
behavior. 

In order to benchmark the four-body energies, Figure 17 compares CG results with fixed-node diffusion Monte 
Carlo (FN-DMC) results carried out by Blume. [i30| From the ground-state energy, the energy crossover curve A4 , 
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FIG. 15: Energy spectrum for four particles with J = in the crossover region (lowest 20 diabatic states). The black curve 
corresponds to the ground state. The blue curves are the states that go diabatically to excited dimer-dimer configurations. 
The red curves correspond to states that connect diabatically to configurations of a dimer plus two free atoms, and the green 
curves correspond to states that connect diabatically to configurations of four free atoms. The lowest green curve is the atomic 
ground state on the BEC side of the resonance. Results from Ref. |129j . 




is constructed as in Refs. [Ml 1130] : 



aw = m^m. (75) 



Here E(A) is the ground-state energy of the four-particle system and E(2) is the ground-state energy of the two- 
particle system. By construction, the energy crossover curve AW varies from 1 in the non-interacting-BCS regime 
(a 0") to in the BEC limit (a — > + ). The energy crossover curve is convenient for comparisons because any 
effec ts o f finite-range interactions on the two-body binding energy are significantly reduced by the subtraction in 
Eq. (75). Therefore, even though both E(A) and E{2) are not completely universal, is universal to a very good 



approximation. 

The solid lines in Figure 17 correspond to the CG prediction while the symbols correspond to the FN-DMC 
predictions obtained by Blume, for an equal mass system. To describe the four-fermion ground state with the FN- 
DMC method, two different guiding wave functions were used: a "normal" and a "paired" guiding function. The 
"normal" wave function can be written as the non-interacting solution multiplied by a Jastrow term that improves 
the description of short-range correlation. This guiding function is well suited to describe the ground state of the 
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system in the BCS-unitarity regime. The "paired" wave function is constructed as an antisymmetrized product of 
two-body solutions and leads to a good description of the ground state in the BEC-unitarity regime. In order to 
reduce finite range corrections, the ranges r of the two-body potentials used in Fig. [17] are set to be much smaller 
than the oscillator lengths, i.e., ro ~ 0.01a? . By analyzing the dependence of the energy on the finite range ro, 
we can extrapolate the zero-range prediction and estimate the deviation in the crossover curve fr om the zero 



range predictions. Such analysis estimates a 1% deviation in the crossover curve presented in Fig. 17 For example, 
at unitarity the CG energies are E — 5.027hu> for r = O.Qlaj^ and E = 5.Q99hu> for r = 0.05a|^T After a more 
careful analysis of the range dependence we obtain an extrapolated zero-range value of E = 5.009huj. For comparison, 
the FN-DMC energy for the square well potential with ro — " 111 
energy calculated by the CG approach. 



O.Ola^ is E = 5.069(9), which agrees well with the 
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FIG. 17: Energy crossover curve A^ as a function of rai^ /a for k = 1. Solid lines are calculated by the CG approach, and 
symbols by the FN-DMC method. Adapted from Ref. [rfb] . 



A test of the validity of the guiding function in FN-DMC calculations is fundamental for an accurate description of 



many-body fermionic systems. Comparisons, such as the one presented in Fig. 17 represent much-needed benchmark 
tests for the "normal" and "paired" trial wave-functions used in the FN-DMC approach |130| The good agreement 

between these two numerical methods suggests that, first, is indeed universal; and second, that both numerical 
methods accurately describe the BCS-BEC crossover for this four-body system. 



F. Extraction of dimer-dimer collisional properties 

In the BEC limit, the lowest four body levels describe different vibrational states of a dimer-dimer configuration. 
Therefore, the systems can be treated effectively in this limit as two-particle systems. A comparison between the two 
particle solutions and the N — 4 solutions allows us to extract information on the effective dimer-dimer interactions. 

To model the effective dimer-dimer interaction we introduce a zero-range pseudopotential. Since the size of the 
dimers are of the order of a, the range of the effective potentials should also be of order of a. Accordingly, effective 
range effects are discussed using an energy-dependent scattering length. Inclusion of the scattering length energy 
dependence is well-known to extend the validity of the zero-range pseudopotential when applied to the scattering of 
two atoms with finite-range potentials under external confinement [1311 1132] . This energy dependence is included 
here through the effective range expansion 



a<id{Ecoi) 
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k 2 r dd . 



(76) 



Here add{E co i) is the energy-dependent dimer-dimer scattering length parameterized by the (zero-energy) scattering 
length add and the effective range r,j_d- The momentum k is associated with the relative kinetic energy of the dimer. 
Thus, k 2 /2fi = E co i where E co i = — 2£ , 2h- The appropriate reduced mass /_t is [idd = M/2, where M is the mass 
of the bosonic molecules, M — mi + rri2- 

Using the effective range expansion [Eq. (65)], the regularized zero-range potential V(r) |133j takes the form 
V{r) — g(E)5(r)(d/dr)r. The scattering strength g is parameterized by the scattering length add and the effective 



33 



range r dd , i.e. 
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The J = spectrum of the two particle trapped system is given by |131] 1132] 
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(78) 



Equation ( 78 ) is a transcendental equation that can be easily solved numerically. The solutions of Eq. ( 78 ) are obtained 



as functions of the a dd and r dd parameters and fitted to the numerical results. The calculation can be carried out at 
different values of the two-body scattering length a and, in this way, one obtains a reliable estimation of a dd and r dd . 

Since the properties of weakly-bound dimers are controlled by the scattering length, it is natural to assume that the 
properties of dimer-dimer interactions are controlled only by a. This implies that a dd and r dd should be proportional 
to the two-body scattering length a. Therefore, we propose expressions of the form a dd = c dd a and r dd — d dd a. The 
parameters c dd and d dd are obtained by fitting the zero-range two-particle solution to the dimer-dimer states of the 
N = 4 system. 
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FIG. 18: Four-body energies of the three energetically lowest-lying dimer-dimer states as a function of a/a^f' for k — 8. 
Panel (a) shows the energetically lowest- lying energy level (i = 0), panel (b) shows the energetically second- lowest (i — 1) and 
panel (c) shows the energetically third- lowest state (i = 2). Circles and crosses show our CG and FN-DMC results, respectively. 
Solid lines show the zero-range model results. Figure from Ref. [130] . 



The results in Fig. 18 illustrate the fitting procedure. The circles in Figs. [I8[a)-(c) show the lowest-lying dimer- 
dimer energy levels, referred to as £?j(4), where i — 0-2 with the center-of-mass energy and the dimer-binding energy 
subtracted. These energies correspond to a two-heavy two-light four-body system with mass ratio k = 8. Solid 
lines represent the energy-dependent zero-range pseudopotential predictions obtained by fitting the two-boson model 
[Eq. ( 78 )] to the four-body energies. The range of atom-atom scattering lengths, a, over which the four-fermion system 
can be described by the two-boson model is significantly extended by the inclusion of the effective range r dd . Such 
inclusion also allows for a more stable and reliable determination of a dd . 

The crosses in Fig. [l8|^a) correspond to FN-DMC energies for the energetically lowest-lying dimer-dimer state. 
The Blume FN-DMC energies, presented in Ref. |130| . are found to be slightly higher than the CG energies, the 
deviation increasing with a. This increasing deviation might be attributed to the variational implications of the 
fixed-node constraint. The functional form of the nodal surface used in the FN-DMC calculations was constructed for 
noninteracting dimers and should be best in the very deep BEC regime. Application of the two-boson model and the 
fitting procedure to the Blume FN-DMC energies provides an alternative determination of a dd and r dd . The increasing 
deviation between the FN-DMC and CG energies with increasing a explains why the effective range predicted by the 
analysis of the FN-DMC energies is somewhat larger than that predicted by the analysis of the CG approach (see 
discussion of Fig. 19 below). 

The formation of few-body states, such as trimers or tetramers, can affect the scattering properties and limit the 
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TABLE II: The dimer-dimer scattering length, add, and dimer-dimer effective range, rdd, obtained using (a) the CG spectrum 
and (b) the FN-DMC energies. The reported uncertainties reflect the uncertainties due to the fitting procedure; e.g., the 
potential limitations of the FN-DMC method to accurately describe the energetically lowest-lying gas-like stateare not included 
here (see Sec. IHB of Ref. [T30]1. 



K 


add/a (a) 


add/ a (b) 


rdd/ a (a) r d d/a (b) 


1 


0.608(2) 


0.64(1) 


0.13(2) 


0.12(4) 


4 


0.77(1) 


0.79(1) 


0.15(1) 


0.23(1) 


8 


0.96(1) 


0.98(1) 


0.28(1) 


0.38(2) 


12 


1.10(1) 


1.08(2) 


0.39(2) 


0.55(2) 


16 


1.20(1) 


1.21(3) 


0.55(2) 


0.60(5) 


20 


1.27(2) 


1.26(5) 


0.68(2) 


0.74(5) 



applicability of the two-boson model. For example, the lowest four-body energy for k = 1, constitutes the true ground 
state of the system, i.e., no energetically lower-lying bound trimer or tetramer states with J n = + symmetry exist. 
However, for larger mass ratios, bound trimer and tetramer states exist. These few-body states can in principle be 
associated with either universal Efimov or nonuniversal physics. Even in the regime where few-body bound states exist, 
the four-body spectrum contains universal states that are separated by approximately 2Huj and are best described 
as two weakly-interacting composite bosons. For fixed a (a > 0), the energy of these "dimer-dimer states" changes 
smoothly as a function of k even in the regime where bound trimer states appear. Thus, our two-boson model and 
the fitting procedure can be extended to this regime of mass ratios. 

Table [TT] and Fig. [19] summarizes the dimer-dimer scattering length and effective range for for selected values of 
the mass ratio K. Circles and crosses in Fig. fl9] correspond to the dimer-dimer scattering length, add, extracted from 
the energies calculated by the CG and the FN-DMC approach, respectively, as a function of k. The add two-boson 
model predictions compare well with those calculated by Petrov et al. within a zero-range framework |134) (solid 



line in Fig. 19). The existence of Efimov states limited Petrov et al. calculations to n < 13.6 since beyond this 



mass ratio a three-body parameter is needed to solve the zero-range four-body equations. As mentioned above, the 
four-body solutions for a finite range potential include deeply-bound solutions which correspond to either a trimer- 
atom or a tetramer. The trimer-atom states affect the dimer-dimer scattering properties only when the dimer-dimer 
configuration is close in energy to the trimer-atom energy, which usually corresponds to a narrow window in two-body 
scattering length. Away from this window, the dimer-dimer states are well described by the two-boson mode which 
predicts a smooth increase of the add up to mass ratio k = 20. 

Note that the existence of energetically-open atom-trimer channels at the dimer-dimer collisional energy provides 
a decay mechanism for the dimers. Clearly, such a decay mechanism cannot be captured within this rather simple 
two-boson model and a more sophisticated treatment should be used to extract the inelastic scattering properties. A 




FIG. 19: Circles and crosses show add /a as a function of k extracted from the four-fermion CG and FN-DMC energies, 
respectively. For comparison, a solid line shows the results from Fig. 3 of Ref. [134] . Diamonds and squares show rdd /a 
extracted from the four-fermion CG and FN-DMC energies, respectively. Figure from Ref. |130| . 



recent theoretical study found good agreement with the two-boson model predictions beyond 13.6 [135] . 

The two-boson model also provides estimates for the dimer-dimer effective range, rdd (shown as diamonds and 
squares in Fig. 19 1. The uncertainty of rdd obtained from the CG approach is estimated to be about 10%; this 
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uncertainty is expected to be larger for the results extracted from the FN-DMC energies since those calculations only 
include one energy curve. As shown in Fig. 19 the ratio r^ja increases from about 0.2 for k = 1 to about 0.5 for 
K = 20. The existence of a finite effective range can be attributed to the effective broad soft-core potential that 
the dimers experience j!34j. Ref. [130] predicts r^d as a function of the mass ratio k. The large value obtained for 
Tdd suggests that effective-range corrections may need to be considered in order to accurately describe the physics of 
molecular Fermi gases. 

To conclude, we have shown that studies of few-body trapped systems can be used to extract information about 
the collisional properties of free systems. Dimer-dimer scattering lengths can be extracted by analyzing the trapped 
few-body spectrum for different two-body scattering length values. Furthermore, energy-dependent corrections to add 
can also be obtained with this method. 



G. Structural properties 



The analysis of the BCS-BEC crossover spectrum can be complemented by an analysis of the wave functions and 
their structural properties. The present section determines the one-body densities and pair-distribution functions for 
two-component Fermi systems in the crossover regime. The averaged radial densities, pi(r), are normalized such that 
47r f pi(r)r 2 dr = 1; Aitr 2 pi(r) represents the probability of finding a particle with mass at a distance r from the 
center of the trap. Here we focus in the mi = m,2 case, where we find that the radial one-body densities pi(r) and 
P2{f) coincide and we can omit the species label. The unequal mass case, in which the radial one-body densities 
Pi(r) and P2(r) differ, will not be considered here but was discussed in Ref. |124j . Also, the averaged radial pair 
distribution functions, Py(r), are normalized so that An J Pij(r)r 2 dr = 1; A^r 1 P.ijir) represents the probability of 
finding a particle of mass rrii and a particle of mass rrij at a distance r from each other. 

These structural properties are computed using the CG method. Since we are only focusing J — states, all the 
structural properties presented here are spherically symmetric. The structural properties are calculated using the 
following general expression: 



Anr 2 F{r) = (V n \8(x - r)|tf„) = / dr 1 ...dr N 8{x - r)\^ n (r u r v , ...,r Nl ,r N2 



(79) 



Here, F(r) is a generic structural property, e.g. the density profiles p± or p2, the interspecies pair-correlation function 
P\2, or the intraspecies pair-correlation functions P\\ or P22- x is the length of the coordinate vector that describes the 
structural property. For p\ and p2, x = n and x = r\>, respectively. For P12, x is the interparticle distance between 
opposite-spin or different species, x — r-yy . For P\\ and P22, x is the same-spin or same-species interparticle distance, 
with x — r\2 and x = ryy, respectively. To evaluate Aitr 2 F(r), ^ n is expanded in the CG basis set, permitting the 
integral in Eq. ( 79 ) to be carried out analytically. 




FIG. 20: Single particle density profiles for interaction strength across the crossover. Thick curves correspond to numerical 
results: a = — aho (solid curves), a — ±00 (dashed curves) and a = O.laho (dash-dotted curves). Thin curves correspond to 
analytic limiting solutions: noninteracting four-fermion prediction (less peaked), dimer-dimer solution in the BEC limit (more 
peaked). 



Figure 



20 



presents the single particle density profiles for different interaction strength across the crossover. The 



3G 



radial density profiles of the equal mass four-fermion system smoothly evolve from the noninteracting solution (less 
peaked thin black curve) in the BCS limit to the prediction of dimer-dimer model in the BEC limit (more peaked thin 
black line). In the BEC limit (a — > + ), the density profile coincides with that of two point like bosonic dimers, i.e., 
p{r) = cxp(— ^V a io )/( a jio v^) 3 - For small but finite positive scattering lengths (dash-dotted curve), the density 
profiles is broader which can be attributed to both the finite size of the dimers and the effective repulsive dimer-dimer 
interaction. 

Next, we analyze the opposite-spin pair distribution function. For a zero-range pseudopotential, the opposite-spin 
pair-correlation function obeys a boundary condition [1361 1137] 



[rPi 2 (r)]U = 2 
[rPi2(r)]r=o a' 

This behavior is a direct consequence of the Bethc-Pcicrls 
[ri 2 4'(ri2)]; i2= o/[ r i2*(' , i2)]r 12 =o = -l/o. The factor of 2 in Eq. 



(80) 

(B-P) boundary condition 
( 80 1 appears because the pair correlation 
function is proportional to the square of the wave function, i.e., Pi2( r ) oc v E'(?'i 2 ) 2 . A direct consequence of Eq. (80) 
is that at unitarity, i.e., when \a\ — oo, r 2 Pi 2 (r) has zero slope at r — > 0. 




FIG. 21: Pair-distribution functions Pia(r), multiplied by r 2 , for equal-mass-two-component Fermi systems with N = 2 
(dashed curve) and N = 4 (solid curve). The N — 2 pair-correlation function has an arbitrary norm selected to match the 
N = 4 pair correlation in the small r regime. 



The numerical verification of Eq. (80) is challenging. For systems with finite range interactions, Eq. (80) is valid in 
a narrow regime of r values. For a < 0, Eq. (801 is valid when r is much larger than the range of the potential and 
much smaller than the mean interparticle distance, i.e., the ro <C r <C a/ lo regime. For a > 0, Eq. (80) is valid when r 
is much larger than the range of the potential and much smaller than the size of the dimer (given by a) and the mean 
interparticle distance, i.e., the ro <C r <C min[a, auo] regime. This regime is almost nonexistent for our numerical 
calculations, because we consider only the a; lo /ro = 100 case. 

The pair distribution function changes considerably as interactions are tuned across the BCS-BEC crossover. Fig- 
ure [22] shows the pair distribution function P\2{r) along the crossover for equal mass systems [k = 1] in the BCS 
(solid lines), unitary (dashed lines) and BEC (dash-dotted lines) regime. As attractive interactions increase, the pair 
distribution function evolves from a single peak structure to a double peak structure. The single peak structure 
is usually associated to solutions where spin-up-spin-down interparticle distances are mainly controlled by a single 
length scale and all particles are roughly at the same distance. This single peak structure appears in the BCS regime 
since the attraction is not strong enough to bind the particles and the typical interparticle distance is set mainly 
by the external trapping potential. The multi peak structures in pair distribution functions are generally associated 
to solutions were more than one length scale control the interspecies interparticle distance. In this case, a double 
peak structure already appears at unitarity and becomes more pronounce in the BEC regime. The peak at small r 
indicates the formation of weakly bound dimers and its width is associated to the size of the dimers. The broader 
peak between 1 a^o and 2 aho is related to the presence of larger atom-atom length scales set approximately by the 
typical dimer-dimer separation in the harmonic oscillator potential. Understanding the BEC solutions as two dimers 
in a trap implies that the four-body configurations includes two interspecies distances of the size of the dimer (short) 
and two interspecies distances of the size of the trap (large). Thus, the probability of finding two particles from 
different species at short distances should be equal to the probability of finding them at large distances. This premise 
can is verified by comparing the area under first and second peaks in the P\2{r)r 2 plot. Finally, note that in moving 
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through the BCS-BEC crossover region, the slope of r 2 P\2(r) at small r changes from positive to zero to negative, as 
is predicted by Eq. (80). 
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FIG. 22: Pair-distribution functions P\i{r), multiplied by r 2 , for equal-mass-two-component Fermi systems with N = 4 and 
J = 0. The solid curve corresponds to a — —dho (BCS regime), dashed curve l/o = (unitarity), and dashed-dotted a — aho 
(BEC regime). Adapted from Ref. 
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FIG. 23: (a) Circles show the pair-distribution function, Pi2(r), multiplied by r 2 , for a = 0.1a; lo (BEC regime). For comparison, 
the blue dash-dotted line shows P\2{r)r 2 for two atoms of mass m with the same scattering length but normalized to 1/2, the 
red dashed line shows Pi2(r)r 2 for two trapped bosonic molecules of mass 2 m interacting through a repulsive effective potential 
with a^d = 0.6a, and the green dotted line shows Pi2(r)r 2 for two trapped noninteracting bosonic molecules of mass 2 m. Panel 
(b) shows a blow-up of the small r region. Figure from Ref. [94] . 



The dimer-dimer model can be quantitatively tested in the deep BEC regime (0 < a <^ a; lo ). In this regime, the two 
peaks of the pair distribution function are well separated and can be independently analyzed. The small-r corresponds 
to dimer formation and is well described by the pair-distribution function multiplied by r 2 and normalized to 1/2, for 
two trapped atoms with a — 0.1a^ o (dash-dotted curve). The large-r peak describes dimer-dimer correlations and is 
well described by the pair distribution function of two bosonic molecules of mass 2m in a harmonic trap (dotted line) . 
The agreement is quite good but it can be improved by including the effective dimer-dimer interaction corrections. 
The dashed curve, almost indistinguishable from the large r part of the pair-distribution function for the four-particle 
system, shows dimer-dimer model prediction including the effective repulsive potential with dimer-dimer scattering 
length add ~ 0.6a |130U134| . Thus, even though the dimer-dimer interaction corrections are small, they are noticeably 
in the pair distribution function. 



1. Dynamics across the BCS-BEC crossover region 



The diabatic representation can be used to ramp an initial configuration through the BCS-BEC crossover region, 
mimicking experiments carried out at different laboratories at JILA and Rice University. The initial configuration is 
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propagated using the time-dependent Schrodinger equation, 



ih d ^=n\{t)m- (si) 

The time dependence of the Hamiltonian comes entirely from that of A(i) = l/a(t) term. In this case, we focus 
on unidirectional ramps. Starting from the ground state on the BCS side, the parameter A is ramped through the 
resonance to the BEC side at different speeds, v = 4r . The relevant dimensionless speed quantity is £ = a,h vfw. The 
parameter £ can be rewritten in terms of is, the density of the system p and the particle mass to relate few-body and 
many-body predictions |138j . This reformulated version of £ (see Ref. 129 ) agrees with the functional form of the 
Landau-Zener parameter obtained in Refs. 139, 140] . The dependence of the Landau-Zener parameter on p has been 
experimentally verified [141j . 

To propagate the initial configuration, we use the diabatic representation obtained previously in Sec |V ET| First, 
we divide the BCS-BEC crossover range into sectors. Starting from the BCS side at A « Xbcs an d finishing at 
the BEC side at A XbeC: the BCS-BEC crossover is divided into 40-80 sectors. At the middle of each sector, 
the time-independent Hamiltonian is diagonalized using the CG method. For four-body systems, thousands of CG 
basis functions are usually needed to describe the spectrum. While in principle this basis set could be used to solve 



Eq. (81 1, in practice this large basis would make the numerical propagation very slow. Instead, we use the diabatic 
representation obtained at the middle of the j-th sector to expand the time-dependent wave function throughout that 
sector, i.e., 



\y(t))=J24(t)W). (82) 

i 

Here, |^) is the diabatic basis function i of sector j, and Nd is the number of diabatic states considered. The time 
dependence only appears in the complex coefficients c-(t). Upon selecting only the lowest 20-100 diabatic states at 



that point, we drastically reduce the size of the Hamiltonian matrix in Eq. (81). Since the inverse scattering length A 



changes very little in each sector, the relevant diabatic states are well described by this reduced basis set throughout 



the sector. The time-dependent Schrodinger equation, Eq. (81 1, is propagated from one edge of the sector to the other 
using an adaptive-step Runge-Kutta method. 

To understand the time propagation of this system, consider the way the probabilities evolve as the system transits 
the BCS-BEC crossover region. At each point of the time propagation, the probability to reside in state i is given by 
Pi(t) = \cl(t)\ 2 . Here j denotes the sector that includes X(t). The probabilities of evolving into a given family can 
be found by summing the probabilities of all states belonging to the same family. For two particles, there are two 
families, the dimcr family, which only includes the lowest state, and the two-atom family which includes the rest of 
the states. In this case, the appropriate definitions are Pd(t) = |c{(i)| 2 , and ffeaC*) = Y^u=2 l°i W| 2 - 

The N = 4 system has four relevant families: the ground state, the excited dimer-dimer states, the dimer-two-atom 
states, and the four-atom states. These families are characterized by the probabilities p g (t), pdd(t), Pd2a(t), and pi a (t), 
respectively. 

Figure [24] presents examples of the numerical time evolution of an N = 4 system during a unidirectional ramp 
at constant speed from the BCS to the BEC side. As expected, the probability of staying in the ground state 
decreases with increasing ramp speed. The probability of evolving into the final four-atom configuration increases 
with speed, which is in agreement with the projection argument. The transfer of probability occurs mainly in the 
strongly interacting regime, —2 < a; lo /a < 2. In this region, the diabatic states are sometimes mixed, producing 
jumps in the probabilities. For example, around a/j Q /a w 1, the red and blue curves have a kink due to an avoided 
crossing between an excited dimer-dimer state and a dimer-two-atom state. 

The probabilities at the end of the time evolution can also be studied as functions of the speed £. Before analyz- 
ing these numerical results, however, consider first the simple Landau-Zener model that provides insights into our 
numerical calculations. 

In its simplest form, the Landau-Zener model, considers a two-level system whose energy difference depends linearly 
on the adiabatic parameter A, i.e., t\ — t2 — &\ and are coupled by £12 independent of A. The time evolution of 
this model can be easily solved, yielding the nonadiabatic transition probability T na to evolve into a final adiabatic 
state different from the initial adiabatic state after the parameter A is varied through an avoided crossing. To obtain 
this probability, the time-dependent Schrodinger equation is propagated starting at t = — oo [with ip(t) = ipi] to 
t = +oo. The nonadiabatic probability is then given by T na = | (^(+oo)|'02) | 2 - Landau and Zener solved this 
problem analytically and showed that 

T na {v) = e~\ (83) 



39 




FIG. 24: Probabilities of different families of the N = 4 system during a unidirectional ramp at constant speed dX/dt from the 
BCS to the BEC side. The initial configuration is cs ). The black curve corresponds to p g {t). The blue curve corresponds to 
Pdd(t)- The red curve corresponds to Pd2a(t), and the green curve corresponds to p4, a {t). (a) Probabilities obtained for ramping 
at a speed of £ w 13. (b) £ w 52. (c) £ « 128 



where 5 = 27ref 2 /(a^). 

Our goal is to use this simple expression to describe each of the important nonadiabatic transitions in the four- 



fermion description. However, equation (83) is not very useful in its current form when it comes to analyzing 
numerical results since it requires a knowledge of nonadiabatic quantities such as a and ei2- £12 can be estimated 
from the difference between adiabatic energy curves at the closest approach. But, in the four-body problem, a is very 
difficult to extract because there is a rich structure of avoided crossings. Nevertheless, Clark [142 showed how a can 
be obtained from an analysis of the P-matrix coupling between two adiabatic states. In the Landau-Zener model, the 
P-matrix between the two adiabatic states \ip+) and \^p~) is 

P + _(A) = (V+fe = z^ TXf vi2 > ( 84 ) 

d\ 4ei2 1 + [a\/{2e 1 2)\ z 

and has a characteristic Lorentzian form whose parameter is directly related to a and ei2- 

Analysis of the spectrum and the P-matrices permits an identification of the states \f/i and involved in the most 
important nonadiabatic transitions, and a determination of the Landau-Zener transition probability 

T ij {v) = e- S * = e~T. (85) 

The Landau-Zener parameter Sij characterizes the transition and is extracted by analyzing the spectrum, while the 
P-matrix is obtained from the numerical description of the adiabatic states using the Hellmann-Feynman theorem. 

Figure [25] (b) displays the results obtained from the numerical time propagation. The black symbols correspond to 
the dimcr-dimcr ground state. The blue symbols correspond to the excited dimer-dimcr family, the red symbols to the 
dimer-two-atom family, and the green symbols to the four-atom family. For slow ramps (small £), the probability of 
forming a "condensate", i.e., remaining in the ground state, is large. For intermediate ramps, the greatest probability 
is to break one bond and end up with a dimer plus two particles. For fast ramps (large £), the probability of staying 
in the atomic ground state on the BEC side is dominant. The probability of the system evolving into an excited 
dimer-dimer configuration remains small for all ramping speeds. 

To analyze these transitions within the Landau-Zener approximation, the partially diabatic states have been labeled 
according to their energies in the BCS regime. This labeling is arbitrary since many of the states are almost degenerate. 
Based on the P-matrix couplings, the possible pathways point towards the states that are most likely important. 
Starting from the ground state, note that has important couplings with states and l^s). Here, is 
the first excited dimer-dimer state, i.e., the lowest state of the excited dimer-dimer family. State ^5) is the first 
excited state of the dimer-two-atom configuration. Since an important probability is transferred to states {^2} an d 
^5), we analyze the couplings of these states to follow the flow of probability. State 1*2) has an important coupling 
with l^s), and state jvE^) has an important coupling with 1^13). Here, ^13) is the lowest state of the four-atom 
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configuration, i.e., the atomic ground state on the BEC side. Figure 25 (a) presents the energy curves of these four 
states. Conveniently, each of these states represents a different configuration. For that reason, this is the minimal 
set of states that can describe the numerical results. While more states could be included in the analysis, here only 
the simplest possible case is considered. The P-matrix analysis also reveals the order in which the transition usually 




FIG. 25: (a) Energy of the most important states amenable to a Landau-Zener description, through the BCS-BEC crossover 
regime. The black curve corresponds to which represents the ground state configuration. The blue curve corresponds to 
|$2) which represents the excited dimer-dimer configuration. The red curve corresponds to l^s) which represents a dimer 
plus two atoms. The green curve corresponds to 1^13} which represents the lowest four-atom configuration, (b) Probability 
of ending up in a given configuration after the ramp, as a function of the dimensionless speed parameter £. The symbols 
correspond to the numerical evolution, while the curves correspond to the Landau-Zener approximation. The colors follow the 
same convention as Figure (a). Results from Ref. |129| . 



occurs. The order of the peaks reveals the following sequence: the first transition is 1 2, then 2 — » 5, then 1—^5, 
and finally 5 13. The Landau-Zener prediction for this sequence is 

Pi = (1-T 1)2 )(1-T 1)5 ), 
P2 = Tx <2 (l - T 2)5 ), 

P5 = ((1 - 7i, 2 )T li5 + T 2i5 T 1)2 )(l - T B) i3), (86) 

P13 = ((1 - 7i )2 )Ti )5 + r 2 ,5Ti.2)r 5 ,i3- 

Again, the sum of all these probabilities is, by construction, unity. The Landau-Zener parameters obtained from the 
P-matrix analysis are 7712 ~ 5.4, 7715 w 6.6, 7725 « 2.1 and 775,13 ~ 13.8. The sequence of Landau-Zener transitions 
in the model shows good agreement with the numerical results even though many possible transitions have been 
neglected in the approximate model. 



2. Universal properties 



The universal properties of two-component Fermi gases interacting through s-wave collisions has been intensively 
studied in recent years. Some particularly important research has been carried out concerning the universal properties 
of four-fermion systems. Here we select two studies that benchmark the universal behavior of four-fcrmion systems. 

The specific point in the strongly-interacting region where the s-wave interaction strength reaches its maximal 
value, is usually called unitarity. Unitarity is alternatively characterized by a divergent s-wave scattering length, 
\a\ = 00. If inelastic two-body scattering channels are energetically open, the scattering length is complex and it does 
not diverge all the way to infinity, but the following discussion assumes that such inelastic processes can be neglected. 
In this situation, if the range of the interaction is much smaller than the typical interparticle distance and if the 
scattering length is divergent, then no relevant length scale exists that can characterize the interaction. This situation 
is similar to the noninteracting limit, where the absence of interactions implies, of course, the absence of a length 
scale that describes the interaction. The absence of a length scale that describes the interaction allows us to extract 
the functional form of various quantities via dimensional analysis. Furthermore, it allows us to relate quantities at 
unitarity to those in the noninteracting limit. 

Dimensional analysis becomes particularly simple in the hyperspherical framework for a noninteracting or unitary 
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system in free space, where the hyperradius is the only coordinate with dimensions of length. Since the potential curves 
have units of energy and the only length scale is given by R, it follows that U(R) oc 1/R 2 . This is equivalent to saying 
that the potential curves at unitarity are proportional to the noninteracting potential curves, i.e.,U(R) oc Uni(R), since 
the noninteracting potential curves have the form 1/R 2 . The resulting predictions have been derived in Refs. |1431I144] . 




FIG. 26: Hyperspherical potential curves at unitarity for the 4-fermion system multiplied by 2/iR 2 /h 2 . The solid lines represent 
the predictions extracted from the spectrum obtained with the CG method. The symbols correspond to direct evaluation of 
the potential curves with the CGHS method. Solid lines correspond to predictions of the large R behavior of the potential 
curves extracted from the analysis of the excitation spectrum of the four-fermion trapped system. 



The four-fermion potential curves calculated at unitarity using CGHS methods allow us to test the premise of 
universality at the four-body level. Figure 26 presents predictions for the lowest 20 adiabatic potential curves for 
a Gaussian interaction model potential. At large R the potential curves become proportional to 1/R 2 as predicted 
by the universal theory. This behavior of the potential curves can be used to understand universal predictions for 
trapped systems such as the virial theorem and the energy spacing of the trapped spectrum. 

The notion of universality extends beyond the unitarity regime and can be applied to any finite scattering length. 
Recently, Tan was able to derive a series of relations between different observables in two-component Fermi sys- 
tems [145H147] . These relations were obtained under the premise of universality and are a consequence of the wave- 
function behavior when two particles come close together (which is well described by the Bethe-Peierls boundary 
condition). The universal Tan relations connect the energy, the expectation value of the trapping potential, the pair 
and momentum distribution functions of two-component Fermi systems through a quantity termed the integrated 
contact intensity 1(a). 
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FIG. 27: Integrated contact intensity obtained in four different ways: through the analysis of the energy dependence on a (solid 
curve), virial theorem (dashed curve), momentum distribution (blue symbols), and pair correlation function (red symbols). 
The solid and dashed curves are essentially indistinguishable on this scale. Figure adapted from Ref. 



Blume and Daily investigated the Tan relations for a trapped four-fermion system [148] . Using the CG method, 
they extracted the spectrum, the pair correlation function, the momentum distribution and the external potential 
expectation value for determining 1(a) in four different ways. Figure 27 presents a comparison of the different predic- 



42 



tions for the contact intensity. The excellent agreement between the different predictions numerically demonstrates 
the validity of the Tan relations. Furthermore, it quantifies the contact intensity for the four-fermion system, which 
can serve as a benchmark for further studies. 



VI. SUMMARY 

This review has concentrated on recent developments in the four-body problem, emphasizing those insights that 
have either used hyperspherical coordinate techniques directly, or else which have benefited indirectly from those 
insights. The utility of formulating the few-body or even the many-body problem in hyperspherical coordinates was 
glimpsed early on by some of the pioneers in the field such as Delves [31| 132]. Smirnov and Shitikova [149j . Macek[38] . 
Lin[3D], Fano |150] . KuppermannjUH], and Aquilanti 55j. This class of methods has been especially valuable for studying 
ultracold collisions in recent years, and ultracold applications have been the focus of this overview. 

Bosonic four-body systems could not be discussed in much detail in this article, owing to length constraints, but 
they have provided an important proving ground for many of the methods discussed in the present review. Early 
ideas on 3-body recombination for three bosons or for three nonidentical fermions that arose first from an adiabatic 
hyperspherical perspective [TUl E] have since been confirmed and developed further, with many useful new insights 
from independent theoretical perspectives. The few-body system receiving the most attention over the years has 
been the three-body problem, both with and without long range Coulombic interactions. The exciting headway 
represented by that body of literature has had some extensions to handle nontrivial reactive processes for systems 
with four interacting particles, 451 and in a handful of studies, for systems with many more particles. }152f[T55] 

Numerous questions still remain to be addressed in this field, aiming in the long run towards not only answering 
questions of universality in the ultracold, but also towards developing systematic ways of handling four-body chemical 
reaction dynamics at the fully quantum level. A large parameter space also remains to be explored just in the ultracold 
limit of the four-body problem, such as varying mass ratios for fermionic and for bosonic systems as well as mixed 
Bose-Fermi systems. 
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Appendix A: Hyperangular coordinates 

As with Jacobi coordinates, there is no unique way to construct the hyperangles of a system. In this section we 
construct the hyperangular coordinates used in the four-fermion problem. The choice of hyperangular parameterization 
has physical meaning. Different parameterizations can be used to describe different correlations within the system. 
Also, in the case of body-fixed coordinates, hyperspherical coordinates can be used to remove the Euler angles of solid 
rotation, reducing the dimensionality of the system. 

1. Delves coordinates 

Unfortunately, or possibly fortuitously depending on your view point, there is no unique way to define the hy- 
perangles in a given system. Here we use a simple, standardized method of defining them used by many others 
[551 |9"TI 11521 1157H159] in the form of the so called Delves coordinates (3TJ [35] . We will begin by examining a well 
known example of hyperspherical coordinates, that of normal spherical polar coordinates. Clearly these coordinates 
can be used to describe the relative motion of 2 particles in 3 dimensions or the position of a single particle in a 
trap-centered coordinate system, but it can also be used in less obvious ways. For instance, spherical polar angles 
may be used to describe the relative motion of 4 particles in 1 dimension. 
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FIG. 28: The tree that gives the standard spherical coordinates for a 3 dimensional system is shown. 




FIG. 29: The canonical tree that gives a hyperangular parameterization for a d dimensional system is shown. 



The components of a three dimensional vector, r, can be written in terms of a radius and two angles as 



x = r cos <p sin t 
y = r sin <fi sin 6 
z = r cos 9. 



(Al) 
(A2) 
(A3) 



This parameterization can be represented in a simple tree structure shown in Fig. |28| The end points of the tree 
represent each component of the vector r, and each node in the tree is represents an angle. Also associated with each 
node is a sub radius. For the lowest node the "subradius" is merely the total length of the vector, r. For the upper 
node the subradius is merely the cylindrical radius p = \/ x 1 + y 2 . Using the tree structure from Fig. 28 a set of rules 
can be developed for extracting the parameterization of Eqs. ( |A1[ ), (A2|, and (A3). Starting at the bottom node with 
total radius, r, move up through the tree to the desired coordinate. For each move through the tree, if you move to 
the left (right) from a node multiply by the sine (cosine) of the angle associated with that node. Continue until you 
reach the Cartesian component. 

This procedure can be generalized readily from three to d dimensions. Start by building a tree with d free ends and 
d — 1 nodes, associate an angle with each node and follow the above rules. Using the tree structure, starting at the 
bottom node with total hyperradius, R, move up through the tree to the desired coordinate. If you move to the left 
(right) from a node, multiply by the sine (cosine) of the angle associated with that node. Continue until you have 
reach the desired Cartesian component. A specific tree for d dimensions is shown in Fig. [29] Following the rules this 
tree gives the hyperangular representation 



d-l 

x n — Rcosa n -i TT sinaj, 

j=n 

< ay < 7T, j = 2,...,d- 1 
< cti < 2tt 



(A4) 
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FIG. 30: The tree structure used to correlate two subspaces to a single hyperradius. 



d-l 

where cosao = 1 an d II sinix,- = 1. This can also be written as 

j=d 



tana„ = — , (A5) 

x n+l 

n = 1,2,3, ...,d- 1. 

This hyperspherical tree has been dubbed the canonical tree [H3 IHZ] as it is simple to construct and very easy to add 
more dimensions to. 

To avoid double counting, the range that the hyperangles take on is restricted depending on how many free branches 
are attached to the node corresponding to a given angle. If the node has two free branches, then the angle takes on 
the full range to 2ir. If the node has one free branch attached, the angle goes from to tt. If the node has no free 
branches attached to it, the associated angle goes from to ir/2. Following these rules for the canonical tree gives 
the ranges of the angles a,, 

< ai < 2tt, 

< ai < n,i — 2, d — 1. 

Another slightly more abstract way of considering this construction is to start by breaking the d dimensional space 
into two subspaces of dimension d\ and cfe, and assuming that these two subspaces are already described by two sets 
of sub-hyperspherical coordinates, {Ri, Qi) and (i?2, fia)- With these assumptions all that remains is to correlate the 
sub-hyperradii. This is done by following the type of procedure described above using the tree structure shown in 
Fig.|3DJ 

i?l =#81110;, (A6) 
i?2 = -Rcosa, 

where a is now the final hyperangle in the system. Using this procedure recursively, one can define the hyperangles in 
the subspaces until the only remaining subspaces are the individual Cartesian components of the total d dimensional 
space. The concept of dividing the total space up into subspaces will prove very or the purpose of constructing basis 
functions. 

As a final example of hyperangular parameterizations, we introduce a parameterization for N 3-dimensional vectors 
{pi} i —i L One could break each vector up into its individual components and use the canonical parameterization from 



Eq. ( A4), but this removes much of the spatial physical intuition that one could bring to bear, such as the individual 
spatial angular momentum corresponding to each vector. Instead one can use a variation on the canonical tree shown 
in Fig. pni On first glance, this tree might seem the same as the canonical tree shown in Fig. [29] In this case, though, 



the large dot at the end of each branch represents the spherical polar sub-tree of the form shown in Fig. 28 for each 
vector pi. Using this tree structure and following the rules outlined above, 2N of the 3A^ — 1 hyperangles are given 
by the normal spherical polar angles for each vector, fa, # 2 , fa, $/v, ^n)- The remaining N — 1 hyperangles are 
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FIG. 31: The tree structure used to parameterize the hyperangles for TV three dimensional vectors is shown. Note that the dot 
at the end of each branch in the tree on the left stands for a spherical polar tree. 



given by 



tan ots 



E - =1 ft 



0<a,< 



Pi+i 

7T 



(A7) 



i = 1,2,3, ...,N- 1. 

where pi is the length of the ith vector. It will be shown in the next section that this parameterization is useful when 
spatial angular momentum plays a role in the problem of interest. For completeness, following the rules laid out in 
Appendix A, the hyperangular volume element that results from this parameterization is given by 



N-i 




cos 2 a-; sin 3 - 7 1 no- 



where 0Ji is the normal spherical polar differential volume for pi. 

The first hyperangular parameterization used here is in the form of Eq. ( A7) for three 3D vectors. The hyperradius 
is defined in the same way as in Eq. (Bl I, 



R 2 



»=i 



(A8) 



where each xi is a Cartesian component of one of the Jacobi vectors. The hyperspherical trees that will be used for 
the four fermion problem will of the type in Fig. [31] in which N vectors are described by the spherical polar angles of 
each vector and a set of hyperangles correlating the lengths of each vector. Specifically, the hyperangles are defined 
by the tree shown in Fig. [32] combined with the spherical polar angles of each Jacobi vector. Following the rules 




FIG. 32: The hyperspherical tree used to parameterize the hyperangular coordinates in the four-fermion problem is shown. See 
Chapter 2 or Ref. [97] for details. 



4G 



described above gives the hyperangles, 



(A9) 



at m , n = tan" 



Here the superscript a — s, il,«2 indicates which Jacobi system is being used, while l,m,n indicate the three Jacobi 
vectors from that system. In principle there are six of these hyperangular coordinate systems that can be constructed 
from each set of Jacobi vectors giving a total of 18 different hyperangular systems. Fortunately we will not need all of 
these to tackle the four fermion problem. If one considers two particles, say 1 and 2, the tree defined by Fig. [32] with 
a = il, I = 2, m = 3 and n = 1 can be decomposed into two subtrees. The right branch describes the hyperangular 
behavior of the dimer alone, while the left branch describes the behavior of the remaining three body system composed 
of a dimer and two free particles. This type of decomposition will be important for evaluating kinetic energy matrix 
elements and defining basis functions. 



Using the recursive definition of the hyperangular momentum operator in Eq. ( B4 ) the total hyperangular momen- 
tum operator can be written in terms of the hyperangular coordinates as 



A 2 =Ai («L,„) 
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(A10) 



where l\ , l m and l n are the normal spatial angular momentum operators for each Jacobi vector. This can also be 
written directly from Eq. (B4) as 



A 2 

A 2 = Ax «„,„) + -tV 
sm ai 
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where all of the hyperangular behavior above the second node in Fig. [32] is described by a sub- hyperangular momen- 
tum, A L,n- _ 

Constructing the hyperspherical harmonics for the four-body system is accomplished following the procedure in 
Section [B] giving 



x N h l j2 sin 1 ' (a; !m )cos' m i a l,m) p ^x + il^-it-l m )/2 ( cos2a i,™) 

X Vlimi 



where P™'^ (x) is a Jacobi polynomial of order 7, y; m (w) is a normal spherical harmonic with spherical polar solid 
angle u, and N^ c is a normalization constant [STll^T] : 



1 y abc 



{2c + d- 



2~j p ^ a+b+c+d+e-2 \ 



c—a—b \ 1 
2 /■ 



1 V 2 



)r 



1/2 



In Eq. (A12| the degeneracy quantum number \i has been replaced with an explicit tabulation of the hyperangular 
momentum quantum numbers, i.e. A/i — > [AA( irn /;, Z m , Z n ]. The total four-body hyperspherical harmonics satisfy the 
eigenvalue equation A 2 Y|^ Zj ; ; , (Q) = A (A + 7) Z; ; l j (fi). The sub- harmonics that are eigenfunctions of 
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Af m can be found as well: 



Y, 



(36) 



X y; itni {uJl)yi m m m (Wm) ■ 



(A13) 



Here the superscript, (3b), indicates that this eigenfunction behaves as a 3-body hyperspherical harmonic. For 
instance if a hyperspherical tree is used with Jacobi vectors defined in the il interaction coordinate system and 1 = 1, 
to = 3, and n = 2, this three-body harmonic describes the free-space behavior of a dimer with two free particles. The 

(n Lm ). The 

(A14) 



three body harmonics obey the eigenvalue equation, A^ m Y"P b ' t[ ; , (f\ m ) = Xi. m (M,m + 4) ^a^'' j, ; 
restrictions on the values of A and A; jm are 

A = Xi, m + l n + 2k 
= h + l m + In + 2j + 2fc, 

where j,k = 0,1,2,.... The quantum numbers li,l m and l n are the spatial angular momentum quantum numbers 
associated with Jacobi vectors pi, p m and p n respectively, and each has a z-projection quantum number associated 
with it which we have suppressed in Eqs. (A12) and (A13). 



,(36) 



2. Democratic coordinates 



Using Delves coordinates greatly simplifies evaluating hyperangular momentum matrix elements, but it still leaves 
the 8 dimensional space of hyperangles. I will only be considering systems with total angular momentum L = 0. 
Therefore it is convenient to move into a body-fixed coordinate system, as the final wavefunction for the four-body 
problem will not depend on the Euler angles that produce a solid rotation of the system. Removing the Euler 
angle dependence is accomplished by transforming into the so-called democratic, or body-fixed coordinates. Four- 
body democratic coordinates are developed in several references (see Refs. [53 [98l |99]). In this work we use the 
parameterization of Aquilanti and Cavalli. For a detailed derivation of the coordinate system see their work in Ref. 



At the heart of democratic coordinates is a rotation from a space fixed frame to a body fixed frame: 

Q = D T (a,/3,7) g bf 



(A15) 



where g is the matrix of Jacobi vectors defined in Eq. (C12), g^f is the set of body fixed Jacobi coordinates, and 
D (a, P, 7) is an Euler rotation matrix defined in the standard way as 



D 



cos a — sin a 
sin a cos a 
1 



cos j3 sin j3 

1 
— sin /3 cos j3 



cos 7 — sin 7 
sin 7 cos 7 
1 



(A16) 



The "~" in Eq. (A15) indicates a transpose has been taken. 



The body-fixed coordinates are defined in a system whose axes are defined by the principle moments of inertia, 
ii, I2 and I3. In this coordinate system the body- fixed Jacobi coordinates are given by 



Qbf 



JXDi 



(A17) 



where D is defined in the same way as in Eq. (A16) with <f>i, </> 2 and ^3 replacing a, ft, and 7. II is a 3 x 3 diagonal 
matrix whose diagonals are given by £1, £2 and £3 which are parameterized by the hyperradius and two hyperangles 
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01 and 02 



* R n 



-R / 

£ 2 = — ^ y 3 sin 2 0! sin 2 2 + cos 2 ©i, 

V *j 

£3 = — ^= \/ 3 sin 2 ©i cos 2 ©2 + cos 2 ©i . 
V3 



(A18) 



To avoid double counting and to allow for different chiralities, ©i and ©2 are restricted to < ©i < ir and < ©2 < 
7r/4. With this parameterization the moments of inertia are given by 



-=C 2 2 +e 3 2 = ^(2 + sin 2 © 1 ) 



A* 
h 



3 — — v 3sin 2 ©! cos 2 2 + 2 cos 2 ©i) , 



(A19) 



R 2 
3 



' 3 = £ 2 + £2 = ~ (3 sin 2 0! sin 2 2 + 2 cos 2 ©1) 



3 



The hyperradius in terms of the principle moments of inertia can be then written as R 2 — £ 2 + £| 
(h + I 2 + h) /2/i. 

With this parameterization, all 8 hyperangles have been defined. The first three are the Euler angles 7}, 
which are external degree s of freedom describing solid rotations of the four body system. The two angles, ©i and 
©2, defined in Eq. (A18|, describe the overall x, y and z extent of the four-body system in the body fixed frame. 
From Eq. (A19), if ©i = 0,7r, then the principle moments of inertia are all equal, i.e. Ji = I2 = I3, meaning that 



the four particles are arranged at the vertices of a regular tetrahedron. When ©i = 7r/2, Eq. (A18) shows that the 
particles are in a planar configuration. The remaining angles, {cf>i, cf>2, ^3}, are kinematic rotations within the system, 
and coalescence points and operations like particle exchange are described in these angles. Broadly speaking, the 
democratic angles ©i and ©2 can be thought of as correlating the overall x, y, and z spatial extent of the four-body 
system in the body-fixed frame, while the kinematic angles <f>i, 4>2, and 03 parameterize the internal configuration of 
the particles. 

Since transformations from one Jacobi set to another are merely rotation matrices (sometimes combined with an 
inversion), the democratic parameterization can always be written in the same form for any given Jacobi coordinate 
system. For example, the symmetry coordinates [Eq. (C3)] can be transformed into the second set interaction 
coordinates [Eq. (Cll)] using the kinematic rotation defined by Eq. (C15). If the democratic parameterization 
defined in Eqs. (A15) and (A17) is used, this transformation reads 



f 2 = D(a,f3,j) nD T i 
= 7) UD T i 



it/. 



>;.2i 



where we have used the fact that the product of two rotations in 3D is itself a rotation. From this it is clear that 
within a given type of Jacobi coordinate (H-type or K-type), all coalescence points arc equally well described. This is 
an important feature as it does not appear in Delves type coordinates, which are strongly dependent on which Jacobi 
tree is used to define them. For the purposes of this paper, to make symmetrization of the wave function easier, we 
will always define the body fixed coordinates in terms of the symmetry Jacobi system q s . 
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Putting this all together, the body-fixed Jacobi vectors in terms of the internal hyperangles can be defined: 

R 

p\ x = —= cos 81 (cos 4>\ cos 4>2 cos 4>3 — sin 0i sin 03) , 
V3 



R I 

ply = —= J sin 2 61 sin 2 6 2 + cos 2 81 (sin 0! cos 2 cos 3 + cos 0i sin 3 ) , 

V o 

I 

p\ z = —= y sin 2 81 cos 2 8 2 + cos 2 81 sin 2 cos 3 , (A20) 

V 3 


p s 2x = —1= cos 81 (cos (f>i cos 2 sin 3 + sin 0i cos 3 ) , 
V3 



P2y = — -p \/ sin 2 81 sin 2 82 + cos 2 81 (sin 0i cos 02 sin 3 — cos 0i cos 3 ) . 
v « 

P2 2 = — 1= ®! COs2 ® 2 + COs2 ®! Sm ^ 2 Sm 03 1 

V d 

/9 3a , = — = cos 81 cos 0i sin 02, 
v3 



P3 y = —= t/ sin 2 61 sin 2 8 2 + cos 2 61 sin 0i sin 2 , 
v 3 



p 3z = — p ^/ sin 2 61 cos 2 62 + cos 2 81 cos 02, 
V 3 

< 9i < vr;0 < 8 2 < j, 

< 01, 02, 03 < 7T. 

The restriction on the range of the internal hyperangles is to avoid double counting configurations and allows for 
configurations of different chirality . 

By moving into democratic coordinates, the dimensionality of the four-body problem can be decreased from 9 to 
6, but, as with many simplifications, there is a cost. This cost comes in the form of the differential volume element 
dn [55]: 

jo (a ■ qjqj \ v / 3cos 3 e2sin 3 e 2 cos2e 2 sin 9 ei 

dil = (dasmpdpdj) pp (A21) 

[(cos 2 8 2 + 3 sin 2 8i cos 2 8 2 ) (cos 2 8 2 + 3 sin 2 8i sin 2 8 2 )] 1 

x <i8id82(i0i sin 02^02^03- 

The first factor is purely from the Euler angle rotation and will always yield a factor of 87T 2 for functions that are 
independent of a, (3 and 7. 

Another price that is paid using democratic coordinates comes in the form of the hyperangular momentum operator, 
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A 2 . In terms of the democratic hyperangles, A 2 is quite complex and can be found in Ref. 

2 (Li + j 2 ; 



A 2 = - A (0i, 8 2 ) + 2fiR 2 ' ll It2 ^ 



2 [h - h) 



I ^- 2 (Ll + Jl) + — — 



h 



2{h-h) 



il-(h-hY 



2 (h - hT 



(Ll + Jl) 



1/2 



+ 



+ 



(h - hY 



1/2 



(h-hY 



-L2J2 



1/2 



(h hY 



L 3 J 3 } , 



where J is the total angular momentum operator, and 
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Terms in J 2 are centrifugal contributions, terms in Li Ji are Coriolis contributions and terms in L 2 are contributions 
from internal kinematic angular momentum with 
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Fortunately, the methods for evaluating matrix elements in what follows will not directly require this form of the 
hyperangular momentum, but it is included here for completeness. 

The final element needed from the democratic coordinates is the inter-particle spacing. The ability to define these 
will be necessary to describe pairwise interactions and correlations. Using Eqs. (CI I, ( C13 1 and (A17), 



1^12 
|»"13 
1^14 
1^23 
1^24 
1^34 



M12 

jr 

Ml4 

K 
M23 

K 

^24 

K 
M23 



11 



22 



st) elf 



22 



22 



(A22) 
(A23) 
(A24) 
(A25) 
(A26) 
(A27) 



■51 



where [ ] f - indicates the ijth element of a matrix. In this equation, only the body fixed Jacobi coordinates from Eq. 
(A17) are used. This is because the unitary Euler rotation used to rotate into the body fixed frame is the same for 
all Jacobi coordinates canceling out the {a,/3,7} dependence in the inter-particle spacings. 

Figure [2] shows the surfaces in {(f>i, (f>2, ^3} for constant in a planar configuration for O2 = 7r/4,7r/6, and n/12 
for equal mass particles. The 4>± coordinate axis has been transformed to cf>i — 7rO (<pi — 7r/2), where (x) is the unit 
step function, to emphasize the symmetry of the surfaces. The red surfaces correspond to the interacting particles 
in the four-fermion system (ri2,ri4,r23 and ^4) while the blue surfaces correspond to the identical fermions (ri3 
and r2i). The identical particle surfaces surround a coalescence point that must be a Pauli exclusion node in the 
final four-body wave function. The simple nature of these coalescence points makes clear the reason for choosing to 
base the democratic coordinates on the symmetry Jacobi vectors. The red surfaces will play an important role in 
the pairwise interaction as these surfaces outline the valleys of the potential. As the system becomes more linear 
(O2 becomes smaller) it can be seen that the surfaces become broader in the </>i direction. In fact, when O2 = (in 
perfectly linear configurations) these surfaces become independent of <f>\. 



Appendix B: Hyperspherical harmonics 

Strictly speaking an overview of hyperspherical coordinates is not really needed for this thesis, as there are many 
excellent existing works on the subject (see for instance [3TJ [3H [!|TJ [S71 11521 1159| V This section is included here for 
completeness in order to provide a more complete foundation for the variational basis functions used in this review. 

To begin, consider a d dimensional Cartesian space whose coordinate axes are given by {xi},,. For the majority 
of this review these coordinates are considered to be the components of a set of Jacobi vectors or the components 
of a set of trap centered vectors, but for now we proceed with a more abstract approach. The basic concept of the 
hyperradius is introduced here as 

d 



r 2 = J2 x I ( B1 ) 



While this definition is used here, often a mass scaling will be inserted. For instance in a trap centered system of 
equal mass atoms an extra factor of 1/N, where N is the number of atoms, will be used to simplify the interpretation 



of the hyperradius. For the purposes of this section, though, this definition will be adequate. With Eq. (Bl|, the d 
dimensional Laplacian can be rewritten in terms of the hyperradius [HI] HZj : 

V 2 = V = 1 9 R 4-i 9 A 2 

^ dx 2 R d ~ l dR OR R 2 ' { 1 

1=1 1 

In this equation A is called the hyperangular momentum, or grand angular momentum operator, the square of which 
is given by 



A 2 = ^-|A y -| 2 , (B3) 

d J)_ 

OXj ax l 

Already the d dimensional Laplacian has a rather pleasing form reminiscent of its 3D counterpart. In fact if d is taken 



to 3, Eq. (51) reduces exactly to the three dimensional Laplacian in spherical coordinates, and A becomes merely the 
normal spatial angular momentum operator. To proceed from here a way of defining the remaining d — 1 degrees of 
freedom in terms of angles is needed. 

The hyperangular momentum operator in terms of hyperangular coordinates can be found by using the fact that 
each subset of Cartesian components is itself a Cartesian vector space. With that in mind, consider the hyperspherical 
tree given by Fig. |30| By writing the Laplacian for each subspace in terms of the sub-hyperradii R\ and R2 and the 
sub- hyperangular momentum operators Ai and A 2 the total hyperangular momentum operator can be extracted |91j . 
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It is 



A 



-1 



,(<il-l)/2, 



, s (d 2 -i)/2 a da 2 



(B4) 



A? + (di - 1) (di - 3) /4 , A 2 + (d 2 - 1) (d 2 - 3) /4 (d-l)(d-3) + l 



sin 2 a 



where a is defined as in Eq. (A6) and A$ is the sub-hyperangular momentum of the subspace of dimension di. If one 



of the subspaces corresponds to a single Cartesian component then the sub-hyper angu lar momentum for that space 
is zero, i.e. if di = 1 then A? = 0. To find Af (A 2 ), ones needs only to apply Eq. (B4| recursively to each subspace. 
In this way, there is a sub-hyperangular momentum operator associated with each node in any given hyperspherical 
tree. 

It is useful to be able to diagonalize the hyperangular momentum operator. The eigenfunctions of A 2 are detailed 
in several references (See Refs. [381 [97 [ [159] for example), and the method of constructing them is given in Appendix 
A. These functions, (17), are called hyperspherical harmonics. Their eigenvalue equation is 



a 2 f Am (n) = a (a + d - 2) (si) 



(B5) 



where A = 0, 1, 2, ... is the hyperangular momentum quantum number. The index \x enumerates the degeneracy for 
each A and can be thought of as the collection of sub-hyperangular momentum quantum numbers that result from a 
given tree. Hyperspherical harmonics are also constructed as to diagonalize the sub-hyperangular momenta of each 
node in a given hyperspherical tree, e.g. 



A 2 1V ( n ) 



Ai (X 1 +d 1 -2)Y Xfl {SI), 



(B6) 



where Ai = 0, 1, 2, ... is the sub-hyperangular momentum quantum number associated with A 2 . The total hyperangular 
momentum quantum number A is limited by the relation 



A = |Ai| + |A 2 



2n, 



(B7) 



where n is a non-negative integer. The absolute values in this case are there to allow for when either d\ or d 2 are 2. In 
this special case the hyperangular momentum quantum number A; associated with the two dimensional subspace can 
be negative, as with the magnetic quantum number, to, in spherical polar coordinates. Equation (B7) only applies if 
both d\ and d 2 are greater than 1. If, for instance d 2 = 1, then the restriction takes on the form 

A = |Ai| + n. 



The behavior illustrated in Eq. (B6| clearly demonstrates why the parameterization shown in Fig. 31 is useful 



Each three dimensional spherical polar subtree will have a spatial angular momentum and z-projection associated 
with it, e.g. 

l\Y Xll (Sl)=k(k + l)Y Xfl (SI), 

where if is the square of the angular momentum operator for the ith vector. This property allows for addition of angular 
momentum in the normal way, through sums over magnetic quantum numbers and Clebsch-Gordan coefficients. Now 
that the hyperspherical harmonics arc defined, it is useful to examine a simple example applying them. 



Appendix C: Jacobi Coordinate systems and Kinematic rotations 



Here we detail the variety of coordinate systems that are used to describe the four-fermion problem in the adiabatic 
hyperspherical framework and the necessary transformations to describe one set in terms of another. The coordinate 
systems used here are not only needed to describe correlations between particles, they also allow the system to be 
reduced in dimensionality by removing the center of mass motion and moving into a body fixed frame. H-type Jacobi 
coordinates are constructed by considering the separation vector for two two-body subsystems, and the separation 
vector between the centers of mass of those two subsystems. K-type Jacobi coordinates are constructed in an iterative 
way by first constructing a three body coordinate set, and then taking the separation vector between the fourth 
particle and the center of mass of the three particle sub-system. 
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1. Jacobi coordinates: H-trees vs K-trees: fragmentation and symmetry considerations 



The first and most obvious symmetry in the four-body Hamiltonian is that of translational symmetry. By describing 
the system in the center of mass frame, the dimensionality of the system can be reduced from d = 12 to d = 9. This 
is done with the use of Jacobi coordinates. In the interest of brevity, We consider here only those coordinates directly 
relevant to the four-fermion problem. These coordinates may be broken into two sets, H-type and K-type, shown 
schematically in Fig. [TJ 

H-type Jacobi coordinates are constructed by considering the separation vector for two two-body subsystems, and 
the separation vector between the centers of mass of those two subsystems, i.e. 



Pcm 




mi 

m,m.. 



f m 2 - 



(CI) 



m k r k + miri \ 
m k +mi J 
m^rA 



TI3 + m 4 

_ (ny +mj) (m t +m k ) 
mi + m 2 + m 3 + m 4 



Here the superscript a enumerates the 24 different H-type coordinates that may be obtained through particle per- 
mutation, p cm is the position of center of mass of the four-body system, and p is an arbitrary reduced mass for the 
four-body system. The prefactors in each Jacobi vector, which are given in terms of the various reduced masses in 
the problem, are chosen to give the so-called mass scaled Jacobi vectors. The kinetic energy in these coordinates can 
be written as 



E 



h2 n2 fi2 n2 

V = V 

2m, r > 2M 



2/i 



where M is the total mass of the four particles. The reduced mass, p, can be chosen to preserve the differential 
volume element for the full 3D problem, ensuring that d 3 p\d 3 p^d 3 p%d 3 p cm = d 3 r 1 d 3 r 2 d 3 r 3 d 3 r i : 



/' = 



mi + m 2 + m 3 + 777,4 



1/3 



Physically, the H-type coordinates are useful for describing correlations between two particles, for example a two body 
bound state or a symmetry between two particles, or two separate two-body correlations. 

When two particles coalesce [e.g. when = rj in Eq. (CI)], the H-type coordinate system reduces to a three body 



system with two of the four particles acting like a single particle with the combined mass of its constituents: 



P? ° = 0, 




m k +m l J 



Locating these "coalescence points" on the surface of the hypersphere is crucial for accurately describing the inter- 
actions between particles, and the coordinate reduction described above will prove useful for the construction of a 
variational basis set. 

K-type Jacobi coordinates are constructed in an iterative way by first constructing a three body coordinate set, and 
then taking the separation vector between the fourth particle and the center of mass of the three particle sub-system, 
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yielding 




(C2) 



mm 



Again a enumerates the 24 different K-type coordinates that result from particle permutations. Examining Fig. [T] 
shows that K-type Jacobi coordinate systems are useful for describing correlations between three particles within the 
four particle system. In the four fermion system, there are no weakly bound trimer states meaning that K-type Jacobi 
coordinates will not be used here, but the methods described in this report can be easily generalized to include these 
type of states. Unless explicitly stated all Jacobi coordinates from here on will be of the H-type, and for notational 
simplicity, we will drop the H superscripts. 



a. Coalescence points and permutation symmetry 

The proper description of coalescence points is crucial for describing two-body interactions, but they are also 
important for describing points of symmetry. For instance if two identical fermions are on top of one another it is 
known that the wave function must vanish at this point owing to the anti-symmetry of fermionic wave functions. 
Here, we are concerned with four fermions in two different "spin" states. Away from a p-wave resonance, the 
interactions between identical fermions can be neglected for low energy collisions. This means that there are two 
types of coalescence points that must be described; two "symmetry" points, when two fermions of the same type are 
on top of each other, and four "interaction" points, places where two distinguishable fermions interact via an s-wave 
potential. 

It might be tempting at this point to choose a single Jacobi coordinate system and then try to describe the 
interactions and symmetries in the same coordinates, but this leads to problems. For instance if it is assumed that 
particles 1 and 3 are spin up and particles 2 and 4 are spin down one might start with coordinates that are simple to 
anti-symmetrize the system in: 



(C3) 

r-2 + r 4 \ 



where it has been assumed that all of the particle masses are equal, mi = m,2 — to 3 = = to leaving \i = m/4 1 / 3 . 
The generalization to distinguishable fermions of different masses is clear. We refer to this Jacobi coordinate system 
as the symmetry coordinates for fairly obvious reasons. If a permutation of two identical fermions is considered, for 
instance 1 and 3, the transformation is simple: 

PuPl = -Pi, (C4) 
PisPl = P2, (C5) 
Pi3Pl = Pi (C6) 
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Similarly for the exchange of particle 2 and 4, 

P2iPl=pl (C7) 

^24P2 = -P2, (C8) 
PvlP% = Pi (C9) 

The points where two identical fermions coalesce are also simply described by taking either p\ — > or p| — > 0. 
The symmetry coordinates are not suited for describing an interaction between two distinguishable fermions, for 
instance 1 and 2. This interaction occurs around the point r x = r 2 . In the symmetry coordinates this means that 



Pi = P2 



V2p s 3 . 



This equation describes a 6 dimensional sheet in the 9 dimensional space, something that is not easy to describe 
directly in any basis set. To get around this problem we introduce two more Jacobi coordinate systems that are useful 
for describing interactions, 




and 



Pi 1 



Pf = \\ — (r a -r 4 ), (CIO) 
; , n+r 2 r 3 + r 4 

P.; 



Pf 



P2 =\Hr to - r 2 ) , (Cll) 



pf 




The superscript il and i2 in Eqs. (CIO I and (Cll) indicate that these Jacobi coordinates are appropriate for inter- 
actions between distinguishable fermions. For instance, a coalescence point between particles 1 and 2 is described 
by P\ — > 0. Another benefit of these coordinates is that they are well suited to describing a dimer wavefunction. If 
particles 2 and 3 are in a weakly bound molecule then the wavefunction for that molecule is only a function of p\ 2 ■ 

Using combinations of these three coordinate systems, pj, p* 1 and p] 2 , can describe all of the possible two-body 
correlations of the fermionic system. This assumes that the system in question is that of four equal mass fermions in 
two internal states with s-wave interactions only. However, the method used is quite general. For instance, K-type 
coordinates can be chosen to describe the possible three-body correlations that can arise due to Efimov physics in 
bosonic systems P21 1160[ 1161 j . 

2. Kinematic rotations 

Since we use different Jacobi systems to describe different types of correlations, a method of transforming between 
different sets of coordinates is needed. In the above section, equal mass particles are considered; extension to arbitrary 
masses is fairly straightforward. To describe the kinematic rotations we keep the masses arbitrary and specify for 
equal masses later. It is convenient here to deal with transforming all of the Jacobi coordinates at once. Thus the 
matrices whose columns are made of the Jacobi vectors are used: 

Q s = {p s nP2,Pl} 

^ = {pf,p 2 \pf} (C12) 

e^ = {p?,p?, P f}. 
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The transformation that takes one coordinate system to another cannot stretch or shrink the differential volume 
element, and thus it must be a unitary transformation. Further, the transformation cannot mix the Cartesian 
components of the Jacobi vector, i.e. pj. has no part of p s y in it. This means that the transformation will be a unitary 
matrix that acts from the right, e.g. 



Q S U S 



(C13) 



The matrices that perform these operations are called kinematic rotations [53 EH [99] , and they will be put to 
extensive use in the calculations that follow. In truth, transformations between coordinates systems that do not 
require an inversion should be considered, but the general principle still holds if improper rotations are included. 
Note that all of the matrix elements must be real, so that the inverse transformation is given merely by the transpose. 

We employ a direct "brute force" method of finding these matrices where the system of equations given in Eq. (CI ) 



are solved for r-y, r 2 , r 3 and in a given Jacobi system. These normal lab-fixed coordinates can then be inserted 
into the definition of the Jacobi coordinates that we wish to describe. The kinematic rotation can then be extracted 
from the resulting relations. Following this procedure gives 
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■'2j 



(C14) 



(C15) 

(C16) 
(C17) 



The same method can be used to find the kinematic rotations to other Jacobi systems, for instance to K-type 
coordinates. 



Appendix D: Implementation of Correlated Gaussian basis set expansion 
1. Symmetrization of the basis functions and evaluation of the matrix elements 

The CG basis functions take the form 

$ A {xi,x 2 , x N ) = S |exp(-ia; T .A.a;)| . 

The symmetrization operator S can be expanded in a set of simple particle permutations, 

N„ 



= Y, sgn(P)|P(A)) 



(Dl) 



(D2) 



Here, N p is the number of permutations that characterize the symmetry 5. Each of these permutations, Pj, has a 
sign associated, sgn(Pj), and is a given rearrangement of the spatial coordinates 



Pi($ A (xi,...,x N ) = $ A (xp i ( 1) ,...,x Pi{N) ) 



(D3) 



The label i characterizes the rearrangement. This rearrangement of the spatial coordinates is equivalent to a rear- 
rangement of the interparticle widths {dij} (or the {ciy}), 



Pfe({<%}) = {dp k Uj}} 



(D4) 
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Therefore, permutation operations can be easily applied and become transformations of the matrix A. 
In general the evaluation of the symmetrized matrix elements of an operator O is, 



(S(A)\0\S(B)) sgn(A)sgn(P^) (P,v(A)|0|p(A)) , (D5) 

i=\ i'=i 

which implies an N p evaluation of unsymmetrized matrix elements. Fortunately, if S(0) = O then, 

(S(A)\0\S(B)) = N p (A\0\S(B)) = N p (S(A)\0\B) . (D6) 



This property significantly reduces the numerical demands since the left hand side of Eq. ( D6 ) implies N p permutations, 
while the right-hand side only implies N p permutations. All operators of the Hamiltoman are invariant under the S, 



operator and their matrix elements obey Eq.(D6|. 

To obtain density profiles and pair-correlation functions, we use the delta function operator. A single delta function 
operator in a given coordinate is not invariant under this transformation; for this reason, the computational evaluation 
is more expensive. Alternatively, we can create a similar operator as a sum of delta functions. If the sum of delta 
functions reflects the symmetry of the problem, the new operator would be invariant under S. 

The permutation operator clearly depends on the symmetry properties of the constituent particles. For identical 
bosons and fermions, 

S = Y^ ai Pi, (D7) 

i=l 

where N p — N\ and a, = 1 for bosons and ai = (— l) p ; p = 0, 1 is the parity of the operator Pj. For two-component 
systems (boson-boson, fermion-fermion, or a Bose- Fermi mixture), 

/v /v 

5=^^« lia ,A4, (D8) 

ii=l i2=l 

where N P1 = Nil, N P2 = A^!, and Ni and N2 are the number of particles in component 1 and 2, respectively. 

The symmctrization operation, if it involves a permutation with a negative sign, can significantly reduce the accuracy 
of matrix elements. In certain cases, the unsymmetrized matrix elements can be almost identical. Because of the 
negative sign of the permutation, the symmetrized matrix elements can become a subtraction of very similar numbers. 
Therefore, accuracy is reduced. These basis functions are usually unphysical, so it is convenient to eliminate them. 
To do this, we evaluate | (S(A)\S(A)) |/ max(| (P i (A)\P i (A)) |). If this is a small number, then the accuracy of the 
matrix elements is reduced. So, in general, we introduce a tolerance of the order of 10 -3 to determine whether to 
keep or discard the basis functions. 

2. Evaluation of unsymmetrized basis functions 

For convenience, we introduce the following simplify notation, 

(cci, • • • , x N \A) = cxp(-^x T .A.x). (D9) 
As a simple example, consider the overlap matrix element 

(A\B) = J dx 1 ..dx N exp{-^x T .{A + B).x). (D10) 
Since the matrix A + B is real and symmetric, there exists a set of eigenvectors y — {yi, ■■.,Vn} with eigenvalues 



/3pf} that diagonalize the matrix. In this set of coordinates, Eq. (D10) takes the simple form, 



(Dll) 
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Here, we used the product pVpV" Av = det(A + B). These basics steps can be followed to evaluate the remaining 
matrix elements. 

To evaluate the kinetic energy, we use the following property, 

{M - \B) = |- (V^AIV.^) . (D12) 

2m 2m 

This property can be simply proven by applying an integration by parts. Also, it simplifies the matrix element 
evaluation and provides an expression which is symmetric in A and B. Then, the matrix element takes the form, 

h 2 N h9 

< A I - 2m ? lB) = 2m 3 Tr((A + ^ ^A.B) (A\B) (D13) 



Another important matrix element, which is similar to Eq.( D13l, is 

(A\x T Cx\B) = ^-3Tr{(A + B^C) (A\B) . (D14) 
2m 

Here, C is an arbitrary matrix. This matrix element is used to calculate the trapping potential energy. In such case, 
C = mu> 2 I /2, where / is the identity matrix. 

Finally, we calculate the matrix element for a two-body central force: 



{A\V{n-r 3 )\B) = J d- i rV(r)(A\5(bi J x-r)\B)=G Clj [V](A\B), (D15) 
where — rj = bfjX, c\- — bfj(A + B)~ 1 b i j, and G C [V] is the Gaussian transform of the potential 

G C [V] = (~) ^ J d\V{v)e-^l\ (D16) 
These matrix elements are enough to describe few-body systems. 



3. Jacobi vectors and CG matrices 



Here, we present the construction of the matrices that characterize the basis functions in terms of the widths 
dij. In the following r = {ri, ...,r^} correspond to Cartesian coordinates, while p = {pi, pjv-i} correspond to 
mass-scaled Jacobi coordinates. First, consider the basis function with the center of mass included 

l-D = * 1 .l/?r.w)-M» | ~E ^^J^ 2 ) =exp(-ir T .A.r). (D17) 

j>i v J 

In the equal-mass case for N particles, it is more convenient to simply use Cartesian coordinates. The ground- 
state ccnter-of-mass wave function of particles in a harmonic trap takes, conveniently, a Gaussian form *&q{Rcm) = 
e -NR 2 CM /2a 2 ho ^ Thus, ^o(Rcm) can be written as ^o(Rcm) = er rT - M ■ r ^ 2 , where M CM is the center-of-mass 
matrix whose matrix elements are M^ M = \/{Na\ ) for all k and Then, for each interparticle distance ry, there 
exists a matrix AfW) so that r?„- = r T .M^.r. The matrix elements of the matrices are m£. = M^P = 1, 

My = Mj % p = — 1; the rest are zero, yielding 

A = M CM + ■ ( D18 ) 

3>i ij 



In some cases it is important to include the center-of-mass motion. For example, this allows one to extract single- 
particle observables such as density profiles. 
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If the center of mass is not included, then Eq. (D18) can be written as, 



cp 



(D19) 



Next, we present the mass-scaled Jacobi vectors (see Fi g. |l| and the corresponding form of the matrices 

Using the H-tree Jacobi coordinates defined in C 1 Eq. <IC lh , the interparticle separation distances can be written: 



r\ - r 2 = pi/di, 



ri - r 3 
r\ - r 4 
r 2 - r 3 
r 2 - r 4 



P3 
P3 
P3 
P3 



Hid 3 

77llG?l 

Hid 3 
midi 
Hid 3 
m 2 di 

Vld 3 



Pi 
Pi 
Pi 
Pi 



m 2 di 
T3 - r 4 = p 2 /d 2 



^ 2 d 3 
m 3 d 2 
H 2 d 3 
mid-2 
H 2 d 3 
m 3 d 2 
H 2 d 3 
m±d 2 



(D20) 
(D21) 

(D22) 

(D23) 

(D24) 
(D25) 



For both the N = 3 and N — 4 systems, the interparticle distances can be written in terms of the Jacobi vectors 

^-^ = E4 ij W (D26) 

k 

Now we can write the matrices M^> in these Jacobi vectors that describe an interparticle distance. The matrix 



elements of these matrices are simply M f 



m 

kl 



4. Selection of the basis set 



There are different strategies for selecting a basis set. If the numbers of dimensions of the system we are studying 
is not that large, then we can try to generate a large basis set that is complete enough to describe several eigenstates 
at different interaction strengths. 

The Gaussian widths dij are selected randomly and cover a range of values from do to the trap length a^ . Specifi- 
cally, the dij are selected randomly using a Gaussian distribution of range 1 and then scaled to three different distances: 
do, an intermediate distance v^oQ/io, and a,ho- These three distances are fixed once the interparticle potential range 
do is fixed. 

The basis set selection depends on the correlation we want to describe. So, the selection process changes depending 
whether the particles are bosons or fcrmions. For fermions, when there is no trimer formation, basis functions with 
more than two particles close together are not important. 

For example, the algorithm for the selection of the basis functions for a two-component four-fermion system divides 
the basis into three parts: the first subbasis generates dij, which are all of the order of atol they are useful for 
describing weakly interacting states. The second subbasis generates two dij of the order of do or y/doaho and the rest 
of the order of a^o! they are useful to describe dimer-dimer states. The third subbasis has one dij of the order of do 
or i/do a /io and the rest of the order of a^ Q . They are useful to describe dimer-two- free- atom states. 



5. Controlling Linear dependence 

A large basis set is usually needed to describe several eigenstates in a wide range of interactions. Since the basis 
set is over complete and the basis functions are chosen semi-randomly, the resulting basis can have linear dependence 
problems. In our implementation we eliminate the linear dependence by reducing the size of basis set. 

To do this, we first diagonalize the overlap matrix and then eliminate the eigenstates with negative or low eigen- 
values. The remaining eigenstates form an orthonormal basis set. Finally, we transform the Hamiltonian to the new 
orthonormal basis set. 



GO 



The threshold for the elimination can be selected automatically taking into account the lowest eigenvalue. If the 
lowest eigenvalue 0\ is small and positive, the tolerance can be selected as, for example, 10 3 Oi. If 0\ is negative 
and the magnitude is large, then the basis set has a lot of linear dependence, and it is more convenient to change the 
initial basis set. 



6. Stochastical variational method 



The SVM has been developed in the context of nuclear physics to solve few-body problems [106 108 . It allows 
a systematically improvement of the basis set. A detailed discussion of the implementation of the SVM will not be 
presented here but can be found in Refs. 1112] , In the following, we present the main concepts of the SVM. 

The SVM is based on the variational nature of the spectrum obtained by a basis set expansion. Consider a basis 
set of size D with eigenvalues {ei,.. ■,££>}, if we add a new basis function then the new eigenvalues {Ai, Ad+i} 
obey Ai < ei < A2...e_o < A^+i. Here, we assume that both set of eigenvalues are arranged in increasing order. 
Thus, by adding a new basis, all the D eigenvalues should decrease or remain the same. Therefore, the lower the new 
eigenvalues are the better the improvement of the basis set. Thus we can test the utility of the added basis function 
by considering the improvement in the eigenvalues. 

In most cases, we are not interested in improving the complete spectrum. To select which states or energies we 
want to improve, we can construct an appropriate minimization function. This function would depend only on the 
energies we want to improve and is minimized by the SVM. 

In order to optimize the basis set, the SVM utilizes a trial an error procedure. Starting from an initial basis set of 
size D, several basis functions are selected stochastically and added, one at a time, to the basis set. For each D + 1 
basis set, the new eigenvalues are evaluated. The basis function that produces the best improvement of the selected 
energies is kept while the remaining basis functions are discarded. The initial basis function is then increased by one 
and the trial an error procedure is repeated. 

If this procedure is continued indefinitely the size of the basis set has become large and the calculations become 
forbiddingly slow. Therefore, it is convenient to increase the basis up to a reasonable size and then continue the 
optimization process without increasing the basis size. This optimization can be carried out by a refinement process. 
Instead of adding a new basis function, we test the importance of the basis functions of the basis set. The trial and 
error procedure is then applied to each of the functions of the basis set. 

For the SVM procedure to be efficient, the evaluation of both the matrix elements and the eigenvalues need to 
be fast and accurate. It is particularly important to obtain very accurate matrix elements because the improvement 
due to a single basis function is usually very small and can only be evaluated reliably if the matrix element are very 
accurate. The matrix element evaluation in the CG and CGHS is both fast and accurate making these methods 
particularly suitable for SVM optimization. 

Also, the evaluation of the eigenvalues can be significantly speeded up in the trial and error procedure. The basis 
functions arc added or replaced one by one which allowing us to reduce the evaluation of the eigenvalues to a root 
finding procedure. This root finding procedure is much faster than any diagonalization procedure. 

The SVM automatically takes care of the selection of the basis function. Also, it tries to avoid linear dependence 
in the basis set by constraining the normalized overlap between any two basis function, i.e., O12/VO11O22, to be 
below some tolerance O max . The tolerance O max is usually selected between 0.95 and 0.99. For example, the size of 
the basis set of N = 3 and 4 can be increased up to 700 and 8000 respectively without introducing significant linear 
dependence. 



Appendix E: Dimer-Dimer Relaxation Rates 



In this appendix I present the derivation of the dimer-dimer relaxation rate used in Section 6.5. This process occurs 
when the two dimers collide causing at least one of the dimers to relax to a deeply bound state. The difference of the 
binding energies is then releases as kinetic energy. This process can be pictured in the hyperspherical picture as an 
infinite series of very closely spaced crossings between the dimer-dimer channel and channels consisting of a deeply 



bound dimer and two free particles. This near continuum of crossings is shown schematically in Fig. 33 'a). 
Using Fermi's golden rule between the initial dimer-dimer state and the final states gives 



v r d i « 1 (*«« (-ft; \y (R, n) I *a (R, n)) | 2 , (ei) 



where ^dd (R', is the dimer-dimer wavefunction, V (R, f2) is the interaction potential and &k is the Ath deeply 
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FIG. 33: (a) A schematic of the channels involved involved inthe dimer-dimer relazaxtion process is shown. The dimer-dimer 
potential (red curve) goes through an inifinite number of crossing with deeply bound states (green dashed curves), (b) The 
hyperradial behavior of the outgoing wavefunction is shown. 



bound dimer state. I now will assume that the dimer-dimer wavefunction is approximated by 

V dd (R; SI) « F dd (R) <5> dd (R; SI) , 



(E2) 



where <& dd (R; SI) is the dimer-dimer hyperangular channel function and F dd (R) is the hyperradial wavefunction 
resulting from the single channel approximation. I further assume that the outgoing deeply bound dimer wavefunction 
can be written as 



(R, n)mip (r 12 ) X (r 3 4, ri 2 ,34) 



(E3) 



where ip (7*12) is the wavefunction for an s-wave deeply bound dimer and Ox (734, ^12, 34) is the free space behavior of 
the resulting three particle system. 



Examining one of the terms from the sum in Eq. ( El ) with a single two-body interaction gives 

2 



V rel 01 



F dd (R) $ dd (R; SI) V 23 (r 23 ) V N X (rk, r 12M ) dRdSl 



(E4) 



where Vff, is the contribution to the relaxation rate by the Ath term in Eq. (El). The first thing to notice is in 



this is that the factor V23 (^23) V 7 ( r i2) is non-zero only when particles 1,2, and 3 are in close proximity, and when 
particles 1, 2, and 3 are in close proximity the remaining degrees of freedom are simplified as well, 



F34 RJ Cri2,34- 

This means that the wavefunction (^34, n.2,34) can be rewritten as 

Ox (r 3 4,f 12 ,34) ~G A (R) f\ (fl) 



(E5) 



(E6) 



where Gx and fx are the hyperradial and hyperangular behavior associated with the Ath outgoing channel. A further 
simplification can be made by realizing that fx (SI) must be independent of SI when particles 1, 2, and 3 are in close 
proximity because the total wavefunction must have zero spatial angular momentum. Rewriting ( E4 1 with these 



G2 



simplifications yields 



Kel « 



/ F dd {R) G k {R) J $ M (R;(l)V 23 (r 23 ) ^ (r 12 ) dQdR 



(E7) 



The hyperangular integral is approximated by the probability that three-particles are close to each other in the 
dimer-dimer channel function. 

And example of G\ is shown in Fig. 33 b) . Away from the classical turning point G\ oscillates very rapidly. This 
fast oscillation will generally cancel out meaning the main contribution to the hyperradial integral is from the region 
near the classical turning point R\. Putting this all together yields 



(E8) 



where T (R\) is the probability that three particles are in close proximity in the dimer-dimer channel function. The 
final step in this derivation is to turn the sum over A into an integral over R\, 



V$oc / p(R x ) 



\F dd (R\ 

rV 



-T(R\)dR) 



(E9) 



where p{R\) is the, nearly constant, density of states. This is possible due to the near-continuum nature of the 
outgoing states. By inserting the WKB approximation wavefunction for F dd (R), the result of Eq. (68) is obtained, 
i.e. 



Pwkb {R\) 
k{R)R 



F(R) dR. 



(E10) 
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